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ABSTRACT 



We present a first-collision model for the evaluation of return flux from the exhaust plume of a 
ring- symmetric HF/DF laser in LEO, generated by an incident flux of ambient molecules traveling at 
orbital speed. The steady plume is bounded by a pair of lip-centered rarefaction fans, and unless 
spacecraft attitude enables incident air molecules to reach the plume through the cavitation regions 
that extend beyond these fans, the spacecraft is shielded from ambient scattering by its own plume. 
Assuming hard-spheres collisions, the first-collision model is given by a simple closed-form expression 
that can be regarded as a source term for scattered exhaust molecules. This source term is integrated 
numerically throughout the fan, yielding the flux arriving at some surface "target point". 
Quantitatively, it is shown that-for a typical HF/DF laser exhaust the contamination level generated 
by ambient scattering is not significant. It was found that the maximum return flux of HF+DF 
constitutes about 2% of the incident ambient flux; this ratio will be nearly constant for LEO 
altitudes. The value of this flux ratio is shown to be dependent on the molecular collision model; it 
may change upon replacing the hard-spheres approximation by a more realistic collision model. A 
possible modification of spacecraft charging by the exhaust was examined, including production of 
HF — and DF ~ . The only significant effect seemed to be shadowing of the downstream half of the 
spacecraft at oblique orbital attitudes. 
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EMPTY PAGE 



1. INTRODUCTION 



This presentation is part of a study on the gas dynamics of ring-symmetric exhaust plumes in 
space, conducted at the Naval Postgraduate School in Monterey. A ring-symmetric jet has zero 
thrust, which makes it suitable as an exhaust configuration for various open loop power plants 
designed to produce high power for relatively short durations. One such system is an envisioned 
space-based chemical laser, shown schematically in Fig. 1-1. In the case of a chemical laser, a ring- 
symmetric configuration would also enable the laser radiation to emerge in the form of an 
axisymmetric beam. 

The exhaust nozzle should be designed to bring the outgoing flow to a supersonic speed at the 
nozzle exit surface. The near field of a free jet is then composed of an inner core bounded by a pair 
of ring-symmetric rarefaction fans centered at the nozzle lips (Fig. 1-1). Beyond the limiting 
characteristic surface of the centered rarefaction waves (CRW), a near-vacuum condition prevails. 
For the purpose of continuum gas dynamic analysis, we assume it is a perfect vacuum. 

An earth orbiting vehicle is subject to an oncoming stream of ambient molecules at a speed of 
U a ~8 (km/ sec), in a direction depending upon its orientation relative to the orbital velocity vector. 
This speed is sufficiently high to cause backscattering of exhaust molecules (see schematic description 
in Fig. 1-2) moving at speeds appropriate to chemical combustion (about 2 to 4 km's). However, 
large exhaust plumes, having achieved stationary flow, may be sufficiently dense at their outer fringes 
to effectively trap and entrain all oncoming ambient molecules. Thus, ambient scattering may be 
significant only in selected ranges of attitude angles, at which ambient molecules can reach the 
vicinity of the spacecraft by traveling almost collisionlessly through cavitation regions. Exhaust 
molecules that may be "candidates" for ambient scattering will hence come from plume segments 
flanked by cavitation regions. The contribution of ambient scattering to contamination will thus be 
highly dependent upon spacecraft geometry and orientation. This may well affect spacecraft design 
and operating procedures. 

The purpose of this report is to present a first-collision model for estimating the flux of exhaust 
molecules backscattered from the fringes of the plume by ambient molecules, along with results of 
sample flux computations performed on a typical HF/DF laser exhaust configuration. The flow field 
throughout the plume is assumed to be governed by the equations of continuum gas dynamics. In 
principle, the liow could be obtained by solving the governing equations, i.e.. the equations for 
stationary isentropic flow in two-dimensional axisymmetric coordinates. In practice, this is normally 
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accomplished by integrating the flow equations in characteristic form, using some finite difference 
scheme (method-of-characteristics). We have performed such computations, but given the complexity 
of applying them to the subsequent integration of ambient scattering flux (due to the need for two- 
dimensional interpolations from an irregular solution grid), we opted for a different alternative : a 
closed-form approximation to the ring- symmetric CRW, based on an analytic expression for flow 
variables along characteristic lines that fan out from the nozzle lip. 

The plan of this report is as follows. In Ch. 2 we outline the approximation to the ring- symmetric 
CRW and present some computation results that demonstrate its accuracy. In Ch. 3 we describe the 
first-collision model and the 3-D spatial integration scheme for computing the flux arriving at the 
cylindrical spacecraft. In Ch. 4 some results of backscattered flux of corrosive molecules (HF + DF), 
showing flux variation with target point location (X g ) and attitude angles (v|/ A , <p A ) are presented. In 
Ch. 5 we take up the subject of spacecraft charging, using results of ambient scattering to assess the 
effect of laser exhaust on spacecraft charging. This is followed by concluding remarks in Ch. 6 and a 
list of references in Ch. 7. A concise description of the flux computation code "AMB" is given in 
Appendix A, followed by the code listing. 
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Figure 1-1. Ring-Symmetric HF DF Laser Exhaust Plume. 




Figure 1-2. Schematic Description of Ambient Scattering. The Cavitation Region is Bounded by 
Lines CA and CB. 
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2. COMPUTATION OF THE PLUME FLOW FIELD 



Most ambient molecules entering the CRW that flanks the exhaust plume are stopped within 
several mean free paths from their point of entry. A quantitative estimate of ambient back-scattering 
would thus depend on the flow field at the outer (hypersonic) fringes of the lip-centered CRW. Even 
though the flow in those regions is generally past the surface of continuum breakdown, the density 
there is reasonably well approximated by the continuum flow field, as demonstrated by Bird's Monte- 
Carlo simulation of a Prandtl-Meyer expansion to vacuum [1] . The evaluation of ambient scattering 
thus calls for an ancillary computational procedure capable of rendering the continuum flow field at a 
large number of points in the ring-symmetric CRW of an exhaust plume. This method was described 
in a recent report [2] . Here we just outline the key ideas and main results of this approximation 
method. 

Our analytic approximation to a ring-symmetric CRW is formulated as follows. In a planar CRW 
(Prandtl-Meyer flow) all flow vanables are uniform along the characteristic lines that fan out from the 
corner (we assume they are the C + family). In the ring- symmetric case the flow near the corner 
approaches asymptotically a corresponding planar CRW flow, which we term the associate CRW. 
However, the gradients along C + characteristics at the corner of a ring-symmetric CRW do not 
vanish as in a planar CRW. The key idea is thus : evaluate flow gradients in C + directions at the 
corner, then use them to extrapolate the associate CRW along C + lines to a finite distance from 
the corner. The extrapolation is a nonlinear function of the radial coordinate y , chosen so that the 
ensuing expression conforms exactly to the flow at the leading (exit) characteristic C + (P,). 
Omitting all details of the analysis, the resulting approximation is presented as the following power- 
law : 

f(«,P) = f(0,P) [y(a,P)/y(0,P)]° ( ° ,P) (2-1) 

w r here f is the streamtube area ratio for isentropic flows ( f= 1 at a sonic point), P is the Mach 
number of a particular characteristic line at the corner, a is a coordinate along the C + (P) 
characteristic line ( a = 0 at the corner), and y is the radial coordinate of a point on the 
characteristic line C + (P) . The Mach number at point (a.P) is readily determined from f(a.P) using 
the standard relation between area ratio and Mach number [3] . A closed-form expression for 5(0. P) 
w'as developed but is not given here; instead, this function is shown in Fig. 2-1. We note that 6 
approaches the asymptotic value of 2/(3-y) as P increases to infinity, and that generally 
1 < 8(0, P) < 2 so that streamtubes diverge at a rate intermediate between that of cylindrical and 
spherical expansion flows. 
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Clearly, in an isentropic flow all thermodynamic variables, and in particular density, can be 
evaluated from f . This approximation is readily applied to the hypersonic portions of a ring- 
symmetric CRW since it turns out that characteristic lines are nearly straight there, which means that 
the characteristic line C + (P) passing through a given point can be readily determined. As a 
demonstration of the degree of accuracy obtainable from this approximation, we show in Fig. 2-2 the 
variation of Mach number along a characteristic line in the ring-symmetric CRW, compared with an 
accurate method-of-characteristics computation. This comparison demonstrates that the analytic 
approximation is reasonably accurate to nearly ten corner-radii away from the comer. 
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3. AMBIENT SCATTERING 



When a rocket or laser exhaust is released into space from an earth- orbiting spacecraft, it 
encounters an oncoming stream of ambient molecules flowing at the orbital speed of U A ^8 (km/sec). 
At altitudes higher than 200 (km), the air/air mean free path exceeds 250 (m), so that it is 
considerably larger than almost any spacecraft. Consequently, ambient molecules would hardly be 
subjected to a shock transition prior to their impact at the spacecraft or exhaust plume. In this 
chapter we describe the formulation of the first-collision model in Section 3.1 and then proceed to 
present the derivation of the flux integration scheme for hard-sphere collisions in Section 3.2. 



3.1 First Collision Model 

The highest ambient number density that we consider for earth-orbiting spacecrafts is 
N A = 1 x 10 16 (m* 3 ), which roughly corresponds to Sunspot Maximum at 200 (km) [4] . The typical 
laser exhaust (Table 4-1) would reach a number density of about 2 x 10 19 (m' 3 ) at the very high Mach 
number of 30. Hence, ambient flux constitutes just a slight perturbation to the near-field portion of 
a typical laser exhaust plume. Obviously, ambient molecules that penetrate the plume, would 
subsequently be entrained by the main flow. But how far do they penetrate? And would exhaust 
molecules scattered by them reach the spacecraft? In seeking answers to these questions, we are led 
to some interesting observations concerning ambient scattering. 

Consider the HF laser depicted in Fig. 1-1. The spacecraft diameter is 5 (m) and the centrally 
located ring-symmetric nozzle is 1 (m) wide. Typical operating conditions (Table 4-1) are assumed. 
They are based on some experimental HF/DF laser studies conducted at TRW [5.6] . Suppose that 
the spacecraft axis is normal to the orbital velocity vector (normal incidence). Let the plane of 
incidence be the plane defined by the intersection of the spacecraft axis with the orbital velocity 
vector. The probability that an ambient molecule traveling in the plane of incidence would reach the 
spacecraft collisionlesslv is exp( — rj), where q is its expected number of collisions with exhaust 
molecules. We define the number q as "molecular thickness", in analogy to "optical thickness". So in 
order to determine the extent to which ambient molecules at normal incidence reach the spacecraft, 
we seek the distribution of radial molecular thickness as function of distance from the spacecraft mid- 
plane (normal to axis at its midpoint). 
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For this purpose we computed the ring-symmetric exhaust flow field, using a semi-inverse 
marching characteristics scheme [7] . The marching was in the radial direction, starting with uniform 
flow at the nozzle exit; the .computation was carried on until it became evident that even at a 
distance of 20 (m) from the mid-plane, the radial molecular thickness was well over 40. The entire 
spacecraft was thus shielded from any ambient scattering at (or near) normal incidence. This 
shielding effect has two significant implications which we discuss briefly below. 

(a) It is present only during stationary’ exhaust flow. At startup and shutdown phases, ambient 
scattering may be substantial even at normal incidence. 

(b) During the stationary phase, ambient scattering is substantial only at attitude angles that enable 
ambient molecules to reach the vicinity of the plume by traveling through "molecularly thin" 
cavitation regions that flank the plume. We thus anticipate a decisive dependence of ambient 
scattering on attitude variations, whenever those variations steer the spacecraft into or out of a 
shielded posture. 

As a first attempt at a quantitative estimate of ambient scattering flux, we have formulated a 
simple first-collision model of this effect. In the sequel we present an outline of the model, along with 
some sample results evaluated for an HF laser configuration identical to that considered for the 
shielding effect mentioned above. 

The basic idea is the following. Ambient molecules entering an exhaust plume, require several 
collisions to become fully "accommodated" with the main flow (i.e., to be entrained by the main flow 
at the prevailing flow velocity and temperature). One may reasonably approximate this process by- 
considering just one collision - the first. 

With the help of some additional assumptions, we were able to derive a closed form expression for 
the flux of exhaust molecules that arrive at the spacecraft following a first collision with an ambient 
molecule. The main assumptions of this model are : 

(1) FIRST COLLISIONS: Only first collisions for either ambient or exhaust molecules are 

considered. Hard-spheres elastic collisions are assumed. Upon a second collision of either an 
ambient or an exhaust molecule, it is considered "lost" (i.e., it joins the main flow). Collisions 
of ambient molecules with spacecraft surfaces are ignored. Ambient molecules are assumed to 
traverse cavitation regions collisionlesslv. 
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(2) COLD FLOW: The oncoming ambient air flow is deemed "cold"; i.e., all molecules move at 

the uniform orbital velocity. The same "cold" assumption is applied to the exhaust flow, since 
most ambient scattering takes place at plume regions of very high Mach numbers (well over 10, 
in the present case). 

(3) CRW Flow Field: ring-symmetric CRW flow field is determined from the power-law 

approximation described in Ch. 2 above. This approximation approaches Prandtl- Meyer flow 
at points whose distance from the nozzle lip is much smaller than the spacecraft radius. 

Based on these assumptions, ambient scattering is represented as a source term for side-scattered 
exhaust molecules, distributed throughout the lip-centered rarefaction fan. The total flux arriving at a 
specified point on the cylindrical spacecraft is readily computed by integrating numerically that source 
distribution over the entire ring-fan. 

The highlights of the spatial integration scheme (Fig. 3-1) are as follows. The limiting 
characteristic surface (M= oo) of the ring-symmetric CRW is divided into surface elements formed by 
dividing the surface into a set of ring-strips which are subdivided in the circumferential (azimuthal) 
direction (<p) into surface elements. The line-of-sight (12) from the "target point" on the spacecraft to 
the center of each surface element is extended into the ring- symmetric CRW, and flux integration 
using the first-collision source term with appropriate weight factors is performed along this line until 
convergence is attained. Contributions from each surface element are summed, taking care to 
disregard portions of the ring-symmetric CRW that are shadowed by the cylindrical spacecraft (either 
the line-of-sight or the trajectory of oncoming ambient molecules may be shadowed). Some further 
details of the flux integration scheme and hard-spheres collisions are provided in Section 3.2 below. 



3.2 Flux Integration Scheme 

The description of the first collision model is hereby supplemented with an outline of the 
expressions used in the flux integration and their derivation. The integration scheme for flux arriving 
at point X s on the spacecraft is depicted in Fig. 3-1. Note that only the plane of incidence is shown 
in Fig. 3-1; at other azimuth angles the geometry is not co-planar, so 3-D geometrical expressions are 
used to get the coordinates (\|/.<p and radial distancete (y 2 + z 2 ) 1 / 2 ) from !2 and S; the derivation of 
these geometrical relations is straightforward, so that we omit these details in the present report. The 
total number flux Q ; (X s ) of i exhaust molecules arriving at point X s is given by the following 
expression : 
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t($) 

%< s > - ? J «ik hi N(n I U(t') - U A | / 1 u A | 



(3-1) 



S 



n ik (S) - ? | dS' <s.. hj N(S0 I u ik (S) - U(S') I / 1 u ik (S) I 



( )l ( )j — Exhaust species 



( ) k - Ambient species 



These expressions are interpreted as follows. The collision depicted in Fig. 3-1 is between exhaust 
molecule mj and ambient molecule m k . The exhaust molar fractions h ; and ambient molar fractions 
h k are assumed uniformly constant, and so are the ambient velocity U A and number density N A - The 
exhaust velocity U(S) and number density N(S) are function of the location in the flow field defined 

<"■* —a 

by £1 and S. These flow variables are computed by first evaluating the coordinates of point QS 
(Fig. 3-1) in the ring- symmetric CRW from the 3-D geometry, and then employing the power-law 
approximation outlined in Ch. 2 above, to get all flow variables for a ring-symmetric CRW. In this 
computation we exploit the fact that characteristic lines fanning out from the nozzle lip are nearly 
straight lines at the low pressure side of the ring-symmetric CRW. 

The £1 integration is performed numerically according to the scheme outlined in Section 3. 1 above, 
as a summation over elements of solid angle (A 3 f2) subtended by area elements on the limiting 
characteristic cone (v|/ = vj/ f ). 

The S integration is considerably more complex. The integrand for this integration is derived as 
follows. Denote by L the line-of-sight distance between point X s and fan point £2.S. A volume 
element at the fan point is given by Av = L" AS A 3 £2. The number of ik pair collisions in Av per unit 
time is <T ;k h. N(S) h k |u(S) — U J exp( — n k (S)| Av , where n k (S) denotes the expected number of 
collisions of ambient molecule k with any exhaust molecule, between its point of entry into the plume 
and point f2,S. We now multiply this term by exp( — r\ ;k (S )| which is the probability that exhaust 
molecule i scattered by ambient molecule k would travel from point £2.S to point X collisionlessly, 
where ri;i<(S) is the expected number of collisions for this path segment. (Note that in Eq. (3-1) the 
summation in the expression for r| ik (S) * s over exhaust species j). 
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The final step in constructing the integrand for the S integration involves the post-collision 

«BaJt 

directional distribution function P ik (S, -£2), whose derivation will be given in the sequel. We 
multiply the integrand by P jk (S,-£2) A £l c which is the fraction of i exhaust molecules scattered by k 
ambient molecules into a solid angle element A 3 £2 C about the unit vector —Q. Considering the flux 
arriving at a surface area element AA $ around point X s , the solid angle element subtended by AA S is 
A 3 n c = AA $ cosa s /L 2 . Eq. (3-1) for Q ; (X s ) now follows upon dividing the resulting expression by 
AA $ , thus referring the arriving flux to a unit area at the point of arrival X s . 

Numerically, the S integration was performed using the classical Runge-Kutta scheme (fourth 
order). The integration for H ik (S) and has to be repeated at each point S. We found 

reasonable convergence with 4 points in the fl k (S) integration and 6 points in the azimuth integration. 
The S integration was terminated when convergence was attained (this is the meaning of the upper 
limit so in the S integral in Eq. (3-1)). The summation over new strips on the limiting cone (\y = y f ) 
was also terminated upon convergence. The CPU time consumed per target point was about 100 
(sec) on IBM 3033 mainframe. 

We now take up the derivation of an expression for the post-collision directional distribution 
function P ;k (S,-£i), which we denote hereafter as P(-£2). We adopt the pair-collision notation 
presented in Fig. 3-2 for the hard-sphere collision analysis. 

As a consequence of conservation of momentum and energy (elastic collisions), the center-of-mass 
velocitv C and the magnitude of the relative velocity C are unchanged bv the collision [8] . The 
post-collision velocities are given by : 

C? = C + p,C* C* = C m ~ ll.C* 

I m r 2 r 2 m "I r 

C r - Cj - c 2 

= mj/On, + m 2 ) 

- O.C, + H,C 2 |c*| = |c,| 

The only free parameter in the expressions for post-collision velocities is the orientation of the 
post-collision relative velocity C . This orientation is uniformly likely to be in any direction in space 
when hard-spheres collision is assumed [8] , as represented by the spherical scattering envelope in 



C r 



C 1 " C ? 



m = m,/(m. +m,) 



(3-2) 
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Fig. 3-3. The probability of obtaining in solid angle element A about — H (Fig. 3-3) is given 
by : 

P(-3) = (i/4k|m 2 cJ 2 ) (aa/a 3 S) = (l/4Jt|cos6|) (|cf| 2 /|p 2 C r | 2 ) (3-3) 

where AA is an area element on the scattering envelope, whose projection on a plane normal to 11 is 
AA |cosS|. We note that the origin of C m in Fig. 3-3 is external to the scattering envelope, resulting 
in two possible scattering elements on the sphere. In all the cases that we computed, however (see 
Ch. 4 below), that point was found to be always internal, so that there was only a single scattering 
solution with post-collision velocity U jk (S) pointing at the spacecraft for any ik pair collision. 
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Figure 3-1. Incidence-Plane Description of Flux Integration Scheme. 
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Figure 3-2. Hard-Spheres Collision Notation. 
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Figure 3-3. Scattering Envelope for Hard-Spheres Collision. 
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4. RESULTS AND DISCUSSION 



We performed several computations of return flux generated by ambient scattering, aimed at 
demonstrating the expected flux level and its variation with spacecraft target point and orbital 
attitude angles. In all these computations we assumed that the exhaust flow is as in the typical 
HF/DF laser case (Table 4-1 below), and that the ambient density and velocity are N A = 1 x 10 16 
(molecules/m 3 ) and U A = 8 (km/sec). As an approximation we further assumed that the sole ambient 
species is molecular nitrogen (molecular weight W A = 28) and that all binary collision cross-sections 
are uniformly given by or=7tD 2 , where D is the molecular diameter (Table 4-1). In each computation 
we evaluated the combined HF + DF flux by assuming that the molar fraction of DF is zero and the 
molar fraction of HF is the combined value for both species (Table 4-1) : .091 + .135 = .226 . This is 
justified by the relatively small difference in molecular weight (just 5%) between these two species. 

Three sets of flux computation were performed as follows : 

(a) Incidence-plane (<p^=0) target points at various distances from the nozzle lip (X s = .l to 
X s = 10 (m)), and at constant incidence angle (\p v = 20°). The results are shown in Fig. 4-1. We 
observe that the flux is fairly insensitive to X s . Also shown in Fig. 4-1 are flux computations 
where the ring-symmetric CRW flow is approximated as a planar CRW (Prandtl- Meyer flow), 
rather than the power-law as in Eq. (2-1) above. The planar case exhibits a somewhat higher 
flux, particularly at large X s . 

(b) Incidence-plane ((p A =0) target points at X s = 1 (m) and at various incidence angles (y A = 0 to 
vp A =40°). A polar representation of the results is given in Fig. 4-2. Note the sharp decrease in 
flux as the incidence angle y A approaches the plume limiting angle \|/ f = 4 1°. 

(c) Azimuth angle variation ( (p ^ = 0 to <p A = 180°) at a constant location (X s =l (m)) and at a 
constant angle of incidence (y 20°). A polar representation of the results is given in Fig. 4-3. 
Observe that flux becomes sensitive to azimuth angle (p ^ only past (p ^ = 90°, where shadowing 
by the cylindrical spacecraft becomes increasingly dominant. 

In addition to return flux we also computed the rms velocity of the arriving molecules. For the 
target points in group (a), the rms velocity varied between 6000 and 6600 (m sec) (the higher velocity 
at smaller X s ), which corresponds to a kinetic energy of about 4 (ev) per molecule (HF). 
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The maximum return flux arriving at the spacecraft is about 0.15 x 10 19 (molecules/m 2 sec), which 
corresponds to a surface deposition rate of about 300 monolayers (HF + DF) per hour. This level of 
contaminating flux may seem to be not outright negligible; however, since return flux is proportional 
to ambient density, it will be scaled down considerably at higher altitudes (and lower ambient 
densities). 

We observe that the maximum return flux constitutes a fraction of about 2% of the incident 
ambient flux. This return flux ratio is roughly maintained at almost all target points and attitude 
angles in groups (a), (b) and (c). The only exceptions are incidence angles near the limiting cone 
(ip = vj/ f ) or at azimuth angles <p A > 125° where shadowing becomes dominant. This observation is 
interpreted as follows. 

Consider the total solid angle subtended by the limiting cone (considered to be infinitely extended 
in the axial direction) as viewed from a target point (for all lines-of-sight £1 pointing outward of the 
cylindrical spacecraft surface). It is independent of target location due to the "self-similar" geometry. 
During each flux computation, we also evaluated the total solid angle subtended by that segment of 
the cone over which the flux integration w'as actually performed (see Section 3.2). It was found out 
that for all but the "shadowed" cases (<p k > 125°), this solid angle constituted a fraction of 86 ±1% 
of the solid angle subtended by the infinite cone. We interpret this result as a hint that geometrical 
"view factors" arising in the course of the flux integration, are not the dominant factor in determining 
the 2% level of flux ratio. What then are the dominant factors? 

For a possible explanation we turn to the flux integration scheme presented in Section 3.2. The 
flux ratio is obtained upon dividing the integrand in Eq. (3-1) by and setting h k = 1 (since we 

assume a single species air). The major factors in the flux ratio integrand appear to be the no- 
collision probabilities exp[ — nj k (S)l and exp[ — fl k (S)j, and the post-collision directional distribution 
function P ;k (S, — £2). The flux-averaged values of these functions in the group (a) computations were 
found to be as follows : P jk (S, — f2) = .09 to .10 , fl ik (S) = .42 to .54 and fl k (S) = .35 to .47. The 
flux-averaged Mach number for group (a) points exhibited a much larger variation : between 30 and 
80, with the higher Mach numbers obtained at further target points. 

These results are interpreted as follows. The ambient no-collision probability exp[ - fl k (S)l is 
sufficiently close to unity, so that in an order-of-magnitude analysis such as the present one, we may 
disregard this factor. If the velocity ratio in the fl ik (S) integral of Eq. (3-1) is assumed to be unity (its 
average value for group (a) points is about 1.4), then the differential in the flux S integration becomes 
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<rN(S)dS — drj ik (S). This implies that the flux S integration results in some average value of the only 
remaining factor : hjP jk (S, — H). Since the H integration introduces a factor of order unity, the 
order-of-magnitude estimate for the arriving-to-incident flux ratio is |hjP jk (S, — £2)] av . The value of 
this estimate is |hjP ik (S, — Q)l av =.226 x .09~.02, which is about equal to the actual flux ratio for 
target points in group (a). 

When an exhaust flow and orbital parameters (velocity and attitude) are specified, P ik (S, — 
depends on the choice of molecular collision model (we chose hard spheres), while hj is uniformly 
constant. The foregoing reasoning thus establishes the collision model as a significant factor in 
determining ambient scattering flux levels, to the extent that P ik (S, — Q) is sensitive to the choice of 
model. 



Table 4-1. Typical Operating Conditions of HF/DF Laser Exhaust 
Mole fractions [H ] = .091 [HF] - .091 [H 2 ] = .104 [DF] = .135 [He] = .579 



Average molecular weight 


7.14 


Specific heats ratio 


1.54 


Stagnation temperature and density 


1400 (K) 

.0075 (kg/m J ) 


Exit Mach number 


4.0 


Molecular diameter (hard spheres) 


2.5x10'^ (m) 


Spacecraft diameter 


5.0 (m) 


Nozzle aperture 


1.0 (m) 



19 




Figure 4-1. Variation of Return Flux with Target Point (X s ). Target Point at Incidence-Plane 
(<P A =0) and Constant Incidence- Angle (y A = 20°). 
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Figure 4-2. Variation of Return Flux with Ambient Incidence Angle (\j/ A ). Fixed Target Point 
(X s = 1 m) Located at Incidence-Plane ((p A = 0). 
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PLANE OF INCIDENCE 



n a = io'W) 




Figure 4-3. Variation of Return Flux with Ambient Azimuth Angle (<p A ). Fixed Target Point 
(X s = 1 m) and Ambient Incidence Angle (\|/ A = 20 o ). 
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5. SPACECRAFT CHARGING 



Spacecraft charging is a major concern to spacecraft designers, particularly for missions in GEO 
and to a lesser extent also in LEO. The exhaust plume of an HF/DF laser operating in the 
ionosphere (300 to 1000 km altitude) may modify significantly the pre-firing charging pattern of the 
spacecraft. Two classes of effects may lead to charging modification; they are: 

(a) The exhaust contains large concentrations of HF and DF molecules which are highly 
electronegative. They may be readily ionized by environmental electrons and change the existing 
spacecraft charging pattern. 

(b) When the spacecraft is oriented obliquely relative to its orbital velocity and the ambient plasma 
impinges at the plume boundary, the plume will cast a "shadow" on the downstream side, 
leading to a very dissimilar charging fluxes on the upstream and downstream halves of the 
spacecraft. 

The knowledge gained in analyzing the ambient scattering effect can be applied to the assessment 
of the effects of ionospheric plasma on spacecraft charging. We first consider the upstream side of 
the spacecraft as mentioned in (a) above. 

We contend that the exhaust-plasma interaction will not drastically alter the charging pattern of 
the upstream half. This assessment is established as follows. Consider the fact that ionospheric 
plasma has a particle number density no higher than 10 12 (m' 3 ) and energy per particle of at most 
1 (ev) (excluding the auroral plasma of polar zones or events of sun storms, where the energy per 
particle is much higher). Significantly, the Debye length at the highest plasma density is very small : 

■y -1 

only about 10' J (m); the largest Debye length in the ionosphere is 10 (m) [9] . Ion thermal velocity 

is typically lower than orbital velocity, but electron velocity is considerably higher than orbital 
velocity (at 1 ev the electron velocity is about U g = 6 * 10 5 m/sec). Hence, ions would typically 
impinge at the plume as a uniform ion beam with the orbital velocity (like ambient molecules), while 
electrons are expected to impinge at the plume with their random-oriented thermal velocity. 

In view of the results of ambient scattering analysis (Ch. 3 and 4 above), and since ions are subject 
to similar collision process with exhaust molecules as neutrals, ions will be stopped at the plume 
fringes much like ambient molecules. By virtue of the small Debye length (typically much smaller 
than the stopping distance), electrons would not penetrate any further than ions, regardless of their 
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collision cross-section with exhaust molecules. The familiar plasma sheath that forms at a solid 
surface, is hence replaced at the plume/plasma boundary by a typically neutral layer whose thickness 
is of the order of an ion/neutral mean free path, but much larger than the Debye length. Only at the 
upper altitude range of the ionosphere does the Debye length become comparable to a plume 
boundary mean free path (about .1 m), but there plasma density and flux are several orders of 
magnitude lower and charging modification is not likely to be significant at the relatively short firing 
duration of about 5 minutes. 

Elastically scattered ions can be deflected towards the spacecraft as a result of elastic collisions 
with exhaust molecules, much like neutrals. Referring to our analysis of the return-to-ambient flux 
ratio (Ch. 4 above), it is clear that the relevant ratio here will be about l/4rt, i.e., of the order of 10% 
(this is due primarily to the role played by the elastic directional distribution function — see Ch. 4). 
A change in the plasma-to-surface current of that order is hence possible, but unlikely to affect 
spacecraft design or operation significantly. The reason is that a design capable of smoothing away 
the inhomogeneous charge flux at oblique attitudes, will not be sensitive to a change in flux pattern 
of the order of 10% (in other words, potential differences may be amplified by 10%, which is hardly 
likely in a sound design to bring about arcing or other threshold phenomena). 

Another effect which may potentially be significant in the upstream half is generation of 
electronegative species (HF — , DF — ) by plasma electrons impinging at the plume. In the sequel, we 
examine the magnitude of this effect, concluding that it is negligible. 

• 

This estimate is best done by considering N , which is the rate of production of HF and DF 

per unit volume, at a typical point in the exhaust where local Mach number is M=30 (this is 
typically the lowest average Mach number for the plume region where ambient scattering takes place 
— see Ch. 4 above). Since energy is released by the electronegative ion formation, the reaction 
involves a third body as follows : 

HF +e — + M -» HF“ + M (5-1) 

where M is the third body molecule. We assume a simplified classical kinetic model for this reaction, 
as follows. The pair HF/M collide with a frequency proportional to the local number density and HF 
molar fraction, and to the average relative velocity. An electronegative ion formation can occur only 
if an electron collides with the pair during their collision, which lasts t c =D/C r , where C r is the 
average relative pair velocity. Based on this classical model, and assuming the same cross-section for 
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electronegative ion formation as for elastic HF/M collisions, the volume rate of electronegative ion 
generation is given by : 

N “ = (tiD 3 N) Nh (JtD 2 U e N e ) (5-2) 

where (JtD 3 N) is the probability that a certain HF or DF molecule will be in contact (D being 

molecular hard-sphere diameter) with any other exhaust molecule (whose number density is N). 

When (7tD 3 N) is multiplied by hN, where h is the HF+DF molar fraction (Table 4-1), the combined 

term reads as the number of colliding HF/M pairs per unit volume. Assuming the electronegative 

formation cross-section is also JtD 2 , the factor 7tD 2 UN where U and N are electron velocitv and 

e e e e 

number density, renders the expression for electronegative generation rate per unit volume. We note 
... • — 

that C r cancels out in deriving Eq. (5-2), so that N does not depend on temperature. This supports 
the use of the kinetic approximation in regions of continuum breakdown (plume fringes are such 
regions). 

• m - 

How is the relative magnitude of N decided? To do that we multiply N by X = l/7tD z N, which 
is the mean free path for a fast moving particle that penetrates the plume. This expression is justified 
by the fact that most incident particles do collide within a distance of order X, and when the particles 

are plasma ions, electrons will adhere to ion spatial distribution by virtue of the small Debye length 

• mmm 

(smaller than X). Thus, XN is the rate of electronegative ion generation per unit area of plume 
boundary'. The ratio P between this rate and the incident electron flux is : 

P“ = XN~/N e U e = (7tD 3 N)h = 2.2 x 10' 10 (5-3) 

where N = 2 * 10* 9 (m‘ 3 ) which corresponds to Mach number M = 30 in the typical case (Table 4-1). 
The fraction of electron flux captured by HF and DF exhaust molecules to form electronegative ions 
is so small (due to the pair-formation term (jtD 3 N)), that it cannot appreciably alter the charging flux 
distribution at the spacecraft surface. 

Another possible effect is the recoil of HF or DF that occurs due to energy released in the 
electronegative formation reaction. The recoiling species might conceivably reach the surface and 
contaminate it. The magnitude of the recoil flux is certainly no larger than P U N = 1300 
(m sec ), where we assume the worst case flux : N e = 10, U e =6 x 10 5 (m/sec) which corresponds 
to about 1 ev energy per electron. This flux level is about 3 x 10* 13 monolayers of HF and DF 
per hour, so that its contribution to surface contamination is utterly negligible. 
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The second kind of charging effects (item (b) above) is due to the fact that the exhaust plume is 
impenetrable to ambient plasma (within a range of sufficiently small distance from the spacecraft, so 
that no extensive diluting of the plume has taken place). The downstream half of the spacecraft in 
oblique attitude will be in the "shadow" with respect to incident plasma. As a first approximation we 
may assume zero plasma flux at the shadowed surface. More accurately, this portion of the 
spacecraft will be subject to a plasma wake flow. However, it is quite difficult to determine the 
charging phenomena that take place in such a wake, as indicated by a recent work on solar sails in 
LEO [9] . Thus, a zero flux at the downstream half seems a practical design assumption. 

Can adverse charging effects occur as a result of shadowing the downstream half? This question 
can be discussed only qualitatively. The reason is that a quantitative analysis requires a lumped- 
circuit model of the spacecraft external surface [10] . Since such a concrete design is not available, we 
can only discuss this question qualitatively. Obviously, assuming zero flux to the downstream half 
during the envisioned 5 minutes of laser firing duration, and requiring that no appreciable voltages 
between the two halves will evolve, leads to the stipulation that the equivalent-circuit 
Capacitance x Resistance should be much smaller than the firing duration. 
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6. CONCLUDING REMARKS 



Our major quantitative conclusion is that for the relatively high ambient density assumed 
(N a = 1 x 10 i6 molecules/m 3 which represents Sunspot Maximum at about 200 km) and for the 
typical HF/DF laser exhaust (Table 4-1), the HF + DF flux backscattered by ambient molecules is 
several hundred monolayers per hour. This flux level may seem as not outright negligible. However, 
since ambient scattering flux is proportional to ambient density, it will be scaled down considerably at 
the lower ambient densities of higher orbital altitudes. 

The operational scenario for HF/DF laser envisions 4 or 5 minutes total operating time; hence the 
contamination by ambient scattering may not be serious due to short operating time. 

The effects of laser exhaust plume on spacecraft charging in the ionosphere were examined. It was 
concluded that the rate of electronegative (HF~ and DF — ) production by impinging electrons was 
negligible; the low rate is a consequence of the assumption that a third body is required to interact 
simultaneously with the HF/e or DF/e pair. No significant modification of charging pattern is 
anticipated. However, at oblique orbital attitudes, the downstream half of the spacecraft will be 
shadowed from the oncoming ambient plasma. This fact has to be reckoned with in designing a ring- 
symmetric laser spacecraft. 

The emphasis in this work was on the method rather than on results. The first-collision model was 
demonstrated to be simple to implement in a code. It is considerably simpler than the more general 
and potentially more accurate Monte Carlo methods commonly used for simulating rarefied flows [8] . 
We found out that the molecular collision model was all important in determining the return flux 
level, which is hardly surprising for scattering by single collision. For the same reason, the collision 
model would also be dominant in a Monte Carlo simulation of the ambient scattering process. 

If and when a mathematical accuracy of the first-collision approximation is established for hard- 
spheres, it might be possible to determine a realistic collision model by comparing computed results 
with measurements. 

This accuracy may be established in either of two ways. One way is by comparison with accurate 
Monte Carlo computations (using hard-spheres collision model). The other is to seek an estimate of 
the error incurred by considering just first collisions and ignoring all subsequent ones. This might be 
achieved by accounting for second collisions in an extended first-collision model, provided a simplified 
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scheme that will obviate the need for increase in the dimensions of the numerical flux integration can 
be devised. We are currently considering such second-collision approaches. 
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APPENDIX A. DESCRIPTION OF AMB CODE 



A.l Description of Subroutines 

We provide a list of the subroutines in the ambient flux integration code AMB for ring-symmetric 

cylidrical spacecrafts. Each subroutine is briefly described. Statements are identified by the 

FORTRAN statement number (columns 1 through 5). 

MAIN PROGRAM The 300 loop is intended to enable several (NCASE) reruns with various data 
in each, all in a single run. Upon calling INI DAT 1, parameters depending on data defined 
in INIDAT are re-computed. The 200 loop is over various XSV(NX) target points. In the 
20 loop the flux integration begins : FLUXC is for particle flux and FLXU2C is for the 
rms of velocity of return flux molecules. All the MAX suffixed parameters denote values at 
which the integrand had the largest value. 

The actual flux integration commences at statement 1 for the summation over strips of 
constant RF. This summation is terminated when convergence is attained (to within 
EPSR). The inner loop 2 is over azimuth angle PHI. Note that the target points are 
generally not in the plane of incidence (PHIA.NE.0), so that no symmetry can be assumed 
in the PHI integration, and it is performed twice in order to cover the entire range in PHI 
(IPAR= l for PHI.GT.0, IPAR = 2 for~PHI.LT.0). The flux integration along the line-of- 
sight is done by calling FLUX. 

INIDAT Initialization of data. There is no input file for this code. INIDAT1 is for parameters 
computed from the data defined by calling INIDAT. 

SOF Stopping routine, called when an error is detected. Here we also trigger a system error by 
computing DSQRT(-l), in order to obtain a calling sequence printout by the operating 
system. 

FLUX This routine calls SUMT for flux integration of one exhaust species at a time. 

LIMIT Here we compute the point of intersection of the line-of-sight with the leading 
characteristic cone. If they do not intersect, the distance of the intersecting point TLIM is 
set to a very large number. 
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SUMT This is the line-of-sight integration routine. Runge-Kutta scheme is used (even though an 
explicit integral is computed). Note that ETAK and ETAIK have to be computed through 
a separate integration at each point of the line-of-sight integration. The integration step 
DT (T=S/RF) is re-adjusted at each integration step. The integration is terminated when 
convergence is attained (to within EPST). 

FETA Here the integrand for the line-of-sight flux integration is evaluated. The hard-spheres 
collision model is used to determine the post-collision directional distribution factor P1K. 
The flux-average of any variable (such as UIK**2 in present version), can be computed by 
summing it multiplied by flux and subsequently dividing by the total arriving flux (see loop 
31 in MAIN PROGRAM). 

PATHIK Here the molecular thickness ETAIK of the I exhaust species scattered by the K ambient 
molecule, is computed by integration along the line-of-sight. 

FT This routine computes the integrand for the ETAIK integration in PATHIK. 

PATHK The analog to PATHIK for K ambient molecule. TAU is the normalized integration 
variable along the trajectory of the penetrating ambient molecule. Note that 
SHADOW = .TRUE, when the trajectory' passes through the cylindrical spacecraft surface 
before entering the fan. 

FTAU Computes the integrand for the ETAK integration in PATHK. 

FAN Computes the fan coordinates PSI, XP, YP for a point on the line-of-sight. It is used to 
determine the Mach number and flow angle from the power-law approximation (see 
MATCH). 

FANT Computes the fan coordinates PSI, XP, YP for a point on the ambient molecule trajectory. 

HMSET Prepares the vector HMV(I) which is the value of the H(M) integral at a set of Mach 
number values (equally spaced in inverse Mach number). This vector is used to compute 
H(M) for an arbitrary M (see HINTER), since this function is needed in the power-law 
approximation of flow in a ring-symmetric fan. Subsequent routines MFUNC, HINTER. 
MATCH and AREAF are all used to implement this approximation. 
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MFUNC Computes the integrand for the H(M) integration in HMSET. 



HINTER Computes H(M) for a given M, from HMV(I) by linear interpolation. Note that the 
interpolation is done with inverse Mach number as the independent variable. 

MATCH Here the approximation to the "inverse problem" of Finding the Mach number at a single 
point in the ring-fan is implemented. An iteration scheme is used to determine the fan 
characteristic passing through the given point [2] . 

AREAF Mach number is computed from value of area ratio function. Newton-Raphson iterations 
are used. 
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A.2 Listing of AMB code 

C$0PTI0NS LIST AMB0001 

C AMBIENT. SCATTERING FROM A RING PLUME BY AMBIENT AIR. AMB0002 

IMPLICIT REAL*8(A-H,0-Z) AMB0003 

COMMON /GAMA/G, Gl, G2, G3, GA, G5, G6, G7, G8, G9, G1 0, Gil, G12,G13,G1 A, G15, AMB 000 A 
1 G16,G17,G18,G19,G20 AMB0005 

COMMON /PAR/CO, ENO, EMI, D, SIGMA, TLIM, DRO , EL 0, QO, TO, FACT, ALOGF, AMB 0006 

1 DPSIO , DTMAX, DETAO , ETALIM, XSI , XSF AMB0007 

COMMON /NPAR/NPHI , IPAR, NP, NR, NX, NXS, NS, NSPEC, NS1 , NS2, NTAUO, NETAO, AMB0008 
1 NAMB, NCASE, ICASE, I FAN AMB0009 

COMMON /GE0M/APF,PAI,PAI2,W,SW,CW,BETA,SBETA,CBETA,PSI1,SPSI1, AMB 0010 

1 CPSI1,PSIF,SPSIF,CPSIF,TPSIF,AK,SK,CK,A0,RF,XF,YF,ZF,AMB0011 

2 PHISOF,PHIF,SPHIF,CPHIF, DYMIN, RMIN,XS,DIST,X0,Y0,Z0, AMB0012 

3 DY0,DEG,PSIN,ST1,CT1, OMEGX, OMEGY , OMEGZ,XSV( 21 ) AMB0013 

COMMON /EPSIL/EPSETA, EPST,EPSR AMB001A 

COMMON /EXT REM/ TEXT (5) , ETA EXT ( 5) , ETAKXT(5) , PHI EXTC5) , AMB0015 

1 PSI EXT ( 5) , EMEXT (5) , FEXT(5) , WEXT(5) , AMB0016 

2 TMAX( 5) , ETAKMX(5) , ETAMAXC5) , PSIMAX(5) , AMB0017 

3 EMMAXC 5) , FMAXC5 ) , AMB0018 

A RFMAXC 5) , PHIFMXC5) , PHIMAXt 5) , WMAX( 5) AMB0019 

COMMON /COUNTS/ICONTC, ICONTT ,ICNT0T,ICNTMX,IQT0T(5),ISHAD(5) AMB0020 

COMMON /SPEC/WAV, XC(5), WC( 5), WRC(5) , XNAME(5) , QFC(5) , QDC( 5), AMB0021 

1 QU2C(5),FLUXC(5), OMEGA( 5) , FLXU2CC 5) , URMSC( 5) AMB0022 

COMMON /AMBI EN/ENA, UA,PSIA,PHIA,HA(3),WA(3), AMB0023 

1 UAX, UAY , UAZ, AA, BA, CA, RA, XA,YA, ZA, SHADOW AMB002A 

COMMON /POINT/XP, YP,XCOR, YCOR AMB0025 

LOGICAL SHADOW AMB0026 

DIMENSION DSUMFC 5) , DSUMDC 5) , DSUMAXC 5) , DSUMU2C5) AMB0027 

NCASE=1 AMB0028 

DO 300 ICASE=1, NCASE AMB0029 

CALL INIDAT AMB0030 

GO TO (301, 302, 303, 30A, 305, 306, 307, 308, 309, 310, AMB0031 

1 311, 312, 313, 31A, 315, 316, 317, 318, 319, 320), ICASE AMB0032 



301 


CONTINUE 


AMB0033 




IFAN=1 


AMB0034 




NXS=3 


AMB0035 




XSI=0 .IDO 


AMB0036 




GO TO 399 


AMB0037 


302 


CONTINUE 


AMB0038 




PHI A=20 . DO/DEG 


AMB0039 




GO TO 399 


AMB0040 


303 


CONTINUE 


AMB0041 




PHIA=50 . DO/DEG 


AMB0042 




GO TO 399 


AMB0043 


304 


CONTINUE 


AMB0044 




PHIA=75. DO/DEG 


AMB0045 




GO TO 399 


AMB0046 


305 


CONTINUE 


AMB0047 




PHI A=1 00 . DO/DEG 


AMB0048 




GO TO 399 


AMB0049 


306 


CONTINUE 


AMB0050 




PHIA=125. DO/DEG 


AMB0051 




GO TO 399 


AMB0052 


307 


CONTINUE 


AMB0053 




PHIA=150. DO/DEG 


AMB0054 




GO TO 399 


AMB0055 


308 


CONTINUE 


AMB0056 




PHIA = 17 5 . DO/DEG 


AMB0057 




GO TO 399 


AMB0058 


309 


CONTINUE 


AMB0059 




GO TO 399 


AMB0060 


310 


CONTINUE 


AMB0061 




GO TO 399 


AMB0062 


311 


CONTINUE 


AMB0063 




GO TO 399 


AMB0064 


312 


CONTINUE 


AMB0065 




GO TO 399 


AMB0066 


313 


CONTINUE 


AMB0067 




GO TO 399 


AMB0068 


314 


CONTINUE 


AMB0069 




GO TO 399 


AMB0070 


315 


CONTINUE 


AMB0071 




GO TO 399 


AMB0072 
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316 


CONTINUE 




AMB0073 




GO TO 399 




AMB0074 


317 


CONTINUE 




AMB 007 5 




GO TO 399 




AMB0076 


318 


CONTINUE 




AMB0077 




GO TO 399 




AMB0073 


319 


CONTINUE 




AMB0079 




GO TO 399 




AMB0080 


320 


CONTINUE 




AMB0081 




GO TO 399 




AMB0082 


399 


CONTINUE 




AMB0083 




PRINT 101 




AMB0084 


101 


FORMAT ( ' 1 ' ) 




AMB0085 








AMB0086 




CALL INDAT1 




AMB0087 








AMB0088 




DO 200 NX=1 , NXS 




AMB0089 




XS=XSV(NX) 




AMB0090 


(X0,Y0,Z0) IS THE POINT AT WHICH FLUX AND DENSITY ARE 


COMPUTED. 


AMB0091 


THE NORMAL TO THE SURFACE AT CX0,Y0,Z0) IS PARALLEL TO 


Y-AXIS. 


AMB0092 




X0 = XS 




AMB0093 




Y0=A0 




AMB 0094 




Z0 = 0. 




AMB0095 




DO 20 NS=NS1,NS2 




AMB0096 




FLUXCC NS ) = 0 . 




AMB0097 




FLXU2C(NS)=0 . 




AMB0098 




OMEGA(NS) = 0 . 




AMB0099 




DSUMAX(NS) = 0 . 




AMB0100 




IQTOT(NS)=0 




AMB0101 




ISHAD(NS)=0. 




AMB0102 




TMAX( NS) =-l . D 44 




AMB0103 




ETAKMXC NS) =-l . D 44 




AMB0104 




PHIMAX(NS)=-1.D 44 




AMB0105 




PHIFMX(NS)=-1 . D 44 




AMB0106 




WMAX( NS) =-l . D 44 




AMB0107 




PSIMAX(NS)=-1.D 44 




AMB0108 




ETAMAX(NS)=-1 . D 44 




AMB0109 




RFMAXC NS) =-l . D 44 




AMB0110 




EMMAXC NS) =-l . D 44 




AMB0111 




FMAX(NS)=-1 . D 44 




AMB0112 


20 


CONTINUE 




AMB0113 




RN=RMIN 




AMB0114 




APF=A0-0.5D0*DR0*SPSIF 




AMB 0 115 




NR= 0 




AMB0116 


1 


NR=NR+1 




AMB0117 




DR=DR0 




AMB0118 




DR=DR0X(APF/A0) 




AMB 01 19 




DR=DR0*(1 ,D0+0.4D0*DR0/XS)**NR 




AMB0120 




RF=RN+DR/2 . DO 




AMB0121 




APF=AO+RF*SPSIF 




AMB0122 




PHISOF=DACOS(AO/APF) 




AMB0123 




DPHI0=0. IDO 




AMB 0 1 24 




NPHI=PHISOF/DPHIO+2 




AMB0125 




DPHI=PHISOF/DBLE(NPHI) 




AMB0126 




DO 21 NS=NS1 , NS2 




AMB0127 




DSUMFC NS ) = 0 . 




AMB0128 




DSUMU2(NS)=0 . 




AMB0129 




DSUMDC NS) =0 . 




AMB0130 


21 


CONTINUE 




AMB0131 




DOMEGR=0 . 




AMB0132 




DO 2 NP=1,NPHI 




AMB0133 




DO 2 IPAR=1,2 




AMB0134 




PHIF=(DBLE(NP)-0.5D0)*DPHI 




AMB0135 




IFCIPAR.EQ.2) PHIF=-PHIF 




AMB0136 








AMB 0137 




CALL FLUX 




AMB0138 








AMB0139 




CROSSl=OMEGY 




AMB0140 




CROSS2=(SPSIF)3«(-OMEGX) + (-CPSIF5(CPHIF)X(-OMEGY) + 




AMB0141 




1 ( -CPSIF^SPHI F)^( -OMEGZ) 




AMB0142 




IFCCROSS1 . LE . 0 . ) 




AMB0143 




1CALL SOF( 'DIRECTION COSINE OF SURFACE NORMAL SHOULD 


BE POSITIVE') 


AMB0144 
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IF(CR0SS2 . LE. 0 . ) 

1CALL SOF( ' NORMAL TO LIMITING CONE HAS NEGATIVE PROJECTION ON LINE 
10F-SIGHT ' ) 

D0MEGA=CR0SS2*DPHI*APF*DR/DIST**2 



DOMEGR=DOMEGR+DOMEGA 
DO 24 NS=NS1,NS2 

DSUMF(NS)=DSUMF(NS)+DOMEGA*QFC(NS)*CROSSl 
DSUMU2(NS)=DSUMU2(NS)+DOMEGA*QU2C(NS)*CROSS1 
IF(DSUMAXCNS) .GT.DOMEGAXQFC(NS)*CROSSl) GO TO 24 
DSUMAX(NS)=DOMEGA*QFC(NS)*CROSSl 
TMAX(NS)=TEXT(NS) 

ETAKMX ( NS )=ETAKXT (NS) 

PHIMAX(NS) =PHIEXT(NS)*DEG 
PHIFMX(NS)=PHIF*DEG 
WMAX(NS)=WEXTCNS)*DEG 
PSIMAX(NS)=PSIEXT(NS)*DEG 
ETAMAX ( NS )=ETAEXT (NS) 

RFMAX(NS)=RF 

EMMAX(NS)=EMEXT(NS) 

FMAX(NS)=QFC(NS)*XC(NS)XQO 
24 CONTINUE 
2 CONTINUE 

DO 26 NS=NS1,NS2 

FLUXCC NS )=FLUXC(NS)+DSUMF( NS ) 

FLXU2CC NS) = FLXU2C( NS) + DSUMU2( NS) 

OMEGAC NS )=OMEGA( NS) + DOMEGR 
26 CONTINUE 
RN=RN+DR 



IF(NR. LE.2) GO TO 1 
IFCNR.GT.99) GO TO 10 
DO 27 NS=NS1,NS2 
IF(FLUXC(NS).EQ.O.) GO TO 27 
ERR=( DSUMF( NS)/FLUXC( NS) )/DOMEGR 
IFCERR.GT.EPSR) GO TO 28 

27 CONTINUE 
GO TO 10 

28 CONTINUE 
GO TO 1 

10 CONTINUE 

DO 31 NS=NS1,NS2 
FLUXC(NS)=XCCNS)*FLUXC(NS)*QO 

OMEGA(NS)=OMEGA(NS)/(2.DOXPAIXDCOS(PSIF/2.DO)X*2) 

FLXU2C(NS)=XC(NS)*FLXU2C(NS)*Q0 

URMSCC NS) =0 . 

IF(FLUXCCNS) .EQ.O.) GO TO 31 
URMSC(NS)=DSQRT(FLXU2C(NS)/FLUXC(NS)) 

C AVERAGE EM (SEE FETA) 

URMSC(NS)= FLXU2C(NS)/FLUXC(NS) 

31 CONTINUE 

PRINT 11,NX,NR,XS,RF, DR, PHISOF*DEG 

11 FORMATC///1X, 'NX,NR,XS,RF,DR,PHIS0F= , ,2I4,3D13.4,F8.4 

1 3X, ' FLUX AND EXTREMA VALUES, ALL SPECIES:’/) 



12 



PRINT 12 

FORMAT (/IX, ’ NAME 

1 ’ FMAX ' 

2 ’ ETAKMX' 

3 ' EMMAX’ 

4 ’ URMSC' 

DO 14 NS=NS1 , NS2 
DLF=0. 



IQTOT ’ , ' ISHAD’, 

’ OMEGA’, ’ TMAX ’ 
’ ETAMAX’,’ PSIMAX’ 
’ RFMAX’ , ’ PI-WMAX’ 
’ FLUXC / LOG'/) 



IF( FLUXCC NS ) . NE . 0 ) 

1DLF=DLOG10(FLUXC(NS))+100 .DO+1 .D-ll 
IDLF=DLF 

DLF=DLF-DBLE( IDLF) 



13 

14 
200 



PRINT 13, XNAMEC NS) ,IQTOT(NS),ISHAD(NS),FMAX(NS), OMEGA ( NS) , 

1 TMAX ( NS) , ETAKMX ( NS ) , ETAMAX ( NS) 

2 PSIMAX ( NS ) , EMMAX ( NS ) , RFMAX ( NS ) 

3 180. D0-WMAX( NS) ,URMSC(NS), 

4 FLUXC(NS) , DLF 

FORMAT (1X,A6,2I6, DIO .3,4F8 .4,4F8 . 1, F8 .2, DIO .3, ’/’,F4.2) 

CONTINUE 

CONTINUE 



AMB0145 

AMB0146 

AMB0147 

AMB0148 

AMB0149 

AMB0150 

AMB0151 

AMB0152 

AMB0153 

AMB0154 

AMB0155 

AMB0156 

AMB0157 

AMB0158 

AMB0159 

AMB0160 

AMB0161 

AMB0162 

AMB0163 

AMB0164 

AMB0165 

AMB0166 

AMB0167 

AMB0168 

AMB0169 

AMB0170 

AMB0171 

AMB0172 

AMB0173 

AMB0174 

AMB0175 

AMB0176 

AMB0177 

AMB0178 

AMB0179 

AMB0180 

AMB0181 

AMB0182 

AMB0183 

AMB0184 

AMB0185 

AMB0186 

AMB0187 

AMB0188 

AMB0189 

AMB0190 

AMB0191 

AMB0192 

AMB0193 

AMB0194 

AMB0195 

AMB0196 

AMB0197 

AMB0198 

AMB0199 

AMB0200 

AMB0201 

AMB0202 

AMB0203 

AMB0204 

AMB0205 

AMB0206 

AMB0207 

AMB0208 

AMB0209 

AMB0210 

AMB0211 

AMB0212 

AMB0213 

AMB0214 

AMB0215 

AMB0216 
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102 

300 



51 



PRINT 102 AMB0217 

FORMATC///1X, 'END RING RUN*,///) AMB0218 

CONTINUE AMB0219 

STOP AMB0220 

END AMB0221 

SUBROUTINE INIDAT AMB0222 

IMPLICIT REAL*8(A-H,0-Z) AMB0223 

REALS LAMDAO , LAMDA1 AMB0224 

CHARACTERS XNAME, XNAMED AMB0225 

COMMON /GAMA/G,G1,G2,G3,G4,G5,G6,G7,G8,G9,G10,G11,G12,G13,G14,G15,AMB0226 

1 G16,G17,G18,G19,G20 AMB0227 

COMMON /PAR/CO, ENO, EMI, D, SIGMA, TLIM,DR0, EL 0,Q0, TO, FACT, ALOGF, AMB0228 

1 DPSIO, DTMAX, DETAO, ETAL IM, XSI , XSF AMB0229 

COMMON /NPAR/NPHI,IPAR,NP,NR,NX,NXS,NS,NSPEC,NS1,NS2,NTAU0,NETA0, AMB0230 

1 NAMB, NCASE, I CASE, I FAN AMB0231 

COMMON /GE0M/APF,PAI,PAI2,W,SW,CW,BETA,SBETA,CBETA,PSI1,SPSI1, AMB0232 

1 CPSI1,PSIF,SPSIF,CPSIF,TPSIF, AK, SK, CK, AO , RF, XF, YF, ZF, AMB0233 

2 PHISOF,PHIF,SPHIF,CPHIF,DYMIN,RMIN,XS,DIST,XO,YO,ZO, AMB0234 

3 DYO,DEG,PSIN,ST1,CT1,OMEGX,OMEGY,OMEGZ,XSV(21) AMB0235 

COMMON /EPSIL/EPSETA,EPST,EPSR AMB0236 

COMMON /EXTREM/TEXT(5),ETAEXT(5),ETAKXT(5),PHIEXT(5), AMB0237 

1 PS I EXT ( 5) , EMEXT ( 5) , FEXT (5) , WEXT ( 5) , AMB0238 

2 TMAX( 5) , ETAKMXC 5) , ETAMAXC 5 ) , PSIMAX ( 5) , AMB0239 

3 EMMAX(5),FMAX<5), AMB0240 

4 RFMAXC 5) , PHI FMX( 5) , PHI MAXI 5) , WMAX( 5) AMB0241 

COMMON /COUNTS/ ICONTC, ICONTT, ICNTOT, ICNTMX, IQT0TC5), ISHAD(5) AMB0242 

COMMON /SPEC/ WAV , XC( 5) , WC(5) , WRC( 5) , XNAMEC5) ,QFC(5),QDC(5), AMB0243 

1 QU2C( 5) , FLUXCC 5) , 0MEGAC5) , FLXU2CC 5) , URMSCC 5) AMB0244 

COMMON /AMBIEN/ENA,UA,PSIA,PHIA,HAC3),WAC3), AMB0245 

1 UAX,UAY,UAZ,AA, BA, CA, RA, XA, YA, ZA , SHADOW AMB0246 

COMMON /POINT/XP,YP,XCOR,YCOR AMB0247 

LOGICAL SHADOW AMB0248 

DIMENSION XCD(5),WCD(5),XNAMED(5) AMB0249 

DATA XCD/.091D0, .091D0, .104D0, .135D0, .579D0/ AMB0250 

DATA WCD/1 .00D0, 20. 0D0, 2. 00D0, 21.000,4.0000/ AMB0251 

DATA XNAMED/’ H • HF H2 ’ DF HE V AMB0252 

DATA IFIRST/O/ AMB0253 

IFAN=2 AMB0254 

PAI=4.D0*DATAN(1.D0) AMB0255 

PAI2=PAI/2 . DO AMB0256 

DEG=180 . DO/PAI AMB0257 

AR=8 . 3143D3 AMB0258 

AV=6 . 022D 26 AMB0259 

OMEGAC=0 . 5 IS FOR HARD SPHERE COLLISIONS, AMB0260 

AN AVERAGE RECOMMENDED VALUE IS ABOUT 0MEGAC=0.75 AMB0261 

OMEGAC=0 . 5D0 AMB0262 

NSPEC=5 AMB0263 

NS1 =2 AMB0264 

NS2=2 AMB0265 

DO 51 NS=1, NSPEC AMB0266 

XCINS) =XCD( NS ) AMB0267 

WC(NS)=WCD(NS) AMB0268 

XNAME(NS)=XNAMED(NS) AMB0269 

CONTINUE AMB0270 

COMBINE HF AND DF MOLE FRACTIONS INTO HF FRACTION AMB0271 

XCC2)=XCC2)+XC(4) AMB0272 

XC( 4) =0 . AMB0273 

AMB0274 

AO = 2 . 5D0 AMB0275 

EM1=4 . DO AMB0276 

RHO0=0 . 0075D0 AMB0277 

T0=1400 . DO AMB0278 

G=1 . 54D0 AMB0279 

D=2.5D-10 AMB0280 

NXS = 1 AMB0281 

XSI=1 . ODO AMB0282 

XSF=10 . DO AMB0283 

AMBIENT AIR AMB0284 

ENA=1 . OOD 16 AMB0285 

UA=8 . D 3 AMB0286 

NAMB=3 AMB0287 

WA( 1 )=28 . DO AMB0288 
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WA( 2) =32 . DO 
WA(3)=16 .DO 
HA( 1)=1 . DO 
HAC2)=0 . 

HA(3)=0 . 

PSIA=20. DO/DEG 
PHIA=O.OODO/DEG 
C INTEGRATION PARAMETERS 
NPHI =6 
NTAUO=4 
NETA0=4 
ICNTMX=100 
RMIN=0 . 

DR0=0 . 1 ODO 
DPSI 0 = 0 . 20D0 
DTMAX=1 .ODO 
DETA0=0 . 50D0 
ETALIM=10 . DO 
EPST=0 . 5D0 
EPSR=0 . 3D0 
FACT=1 . D 20 
RETURN 

CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 
C COMPUTATION OF DATA-DEPENDENT PARAMETERS 

CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX*X 

ENTRY INDAT1 

CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 

ALOGF=DLOG< FACT) 

WAV=0. 

DO 52 NS=1 , NSPEC 
WAV=WAV+XC(NS)XWC(NS) 

52 CONTINUE 

DO 53 NS=1, NSPEC 
WRC( NS) =WC( NS)/WAV 

53 CONTINUE 
SIGMA=PAIXDXX2 
ENO=RHOOXAV/WAV 
CO=DSQRTCGXARXT O/WAV) 

XSV(1)=XSI 

IF(NXS.EQ.l) GO TO 12 

DXL=(DL0G(XSF)-DL0G(XSI))/(DBLE(NXS)-1 .DO) 

XLI=DLOG(XSI ) 

DO 11 NX=2,NXS 

XSV( NX) = DEXP( XLI + ( DBL E(NX)-l.DO)XDXL) 

11 CONTINUE 

12 CONTINUE 

G1 = (G-1 . D0)/2 . DO 
G2=(G+1.D0)/(2.D0X(G-1.D0)) 

G3=G/2 . DO 

G4=(G+1.D0)/(G-1.D0) 

G5=DSQRT((G+1.D0)/(G-1.D0)) 

G6=l .D0/(G-1 .DO) 

G7=2.D0/(G+1 .DO) 

G8 = ( 0 . 5D0x(G+1 . D0)x* 2/(G-1 . DO))xx(l . DO/CG+l . DO) )X 
1 ((G+1.D0)/(G-1.D0))XX((G-1.D0)/(G+1.D0)) 

G9=(G+3.D0)/(2.D0X(G-1 .DO)) 

G10=(7.D0-3.D0XG)/(2.D0X(G-1.D0)) 

G11=(5.D0-3.D0XG)/(2.D0X(G-1 .DO)) 

G13= ( 2 . DO-G)/( 2 . D0X(G-1 . DO) ) 

G14=G/(2.D0X(G-1.D0)) 

G15=(G+l.D0)/(3. DO-G) 

ZETA1=G5XDATANCDSQRT(EM1XX2-1 .D0)/G5) 

AMU1=DASIN(1 .DO/EMI) 

PSI1=PAI2+AMU1 

SPSI1=DSIN(PSI1) 

CPSI1=DC0S(PSI1) 

PSIF=PAI2+AMU1+ZETA1-G5XPAI2 

SPSIF=DSIN(PSIF) 

CPSIF=DCOS(PSIF) 

TPSIF=DTAN(PSIF) 

TETA1=PSI1-AMU1 

ST1=DSIN(TETA1) 



AMB0289 
AMB0290 
AMB0291 
AMB0292 
AMB0293 
AMB0294 
AMB0295 
AMB0296 
AMB0297 
AMB0298 
AMB0299 
AMB0300 
AMB0301 
AMB0302 
AMB0303 
AMB0304 
AMB0305 
AMB0306 
AMB0307 
AMB0308 
AMB0309 
AMB0310 
AMB0311 
AMB0312 
AMB0313 
AMB0314 
AMB0315 
AMB0316 
AMB0317 
AMB0318 
AMB0319 
AMB0320 
AMB0321 
AMB0322 
AMB0323 
AMB0324 
AMB0325 
AMB0326 
AMB0327 
AMB0328 
AMB0329 
AMB0330 
AMB0331 
AMB0332 
AMB0333 
AMB0334 
AMB0335 
AMB0336 
AMB 03 37 
AMB0338 
AMB0339 
AMB0340 
AMB0341 
AMB 0342 
AMB0343 
AMB0344 
AMB0345 
AMB0346 
AMB0347 
AMB0348 
AMB0349 
AMB0350 
AMB0351 
AMB0352 
AMB0353 
AMB0354 
AMB0355 
AMB0356 
AMB0357 
AMB0358 
AMB0359 
AMB0360 
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CT1=DC0S(TETA1) 

QO=ENA*UA 

LAMDA0=1 . D0/(DSQRT(2.D0)*SIGMA*EN0) 
LAMDA1=LAMDA0*( 1 . D0+G1*EM1**2 )**(G6 -OMEGAC+O . 5D0 ) 
AA=DCOS(PSIA) 

BA=DSIN(PSIA)XDCOS(PHIA) 

CA=DSIN(PSIA)*DSIN(PHIA) 

UAX=-UA*AA 
UAY=-UA*BA 
UAZ=-UA*CA 
XCOR=0 . 

YCOR=AO 



NSPEC=’,I3/ 

’,11(2X,A6,2X)) 



PRINT 201 , NSPEC, XNAME 

201 F0RMAT(/1X, 'SPECIES DATA 

1 IX, 'SPECIES NAMES 

PRINT 202 XC 

202 FORMAT ( IX, 'MOLE FRACTION XC= • , 11 ( F8 . A, 2X) ) 

PRINT 203, WC 

203 FORMAT ( IX, ’MOL. WEIGHT WC= • , 11 (F8 . 4, 2X) ) 

PRINT 21,AR,AV,WAV,G,RHOO,TO,ENO,CO,D 

21 FORMATC/1X, ’THERMODYNAMIC DATA'/ 

1 IX, ' AR, AV, WAV , GAMMA = ',2X,2D14.5,2F9.3/ 

2 IX, ’RHOO,TO,ENO,CO,D=',D12.4, F8 . 0, D13 . 5, 2D12 . 4) 

PRINT 22,EM1,PSI1*DEG,PSIF*DEG, 

1 AO , LAMDAO , LAMDA1 

22 FORMATC/1X, 'FLOW AND GEOMETRY DATA’/ 

1 IX, ’EM1,PSI1,PSIF=',3F9.3/ 

2 IX, ’ AO , LAMDAO , LAMDA1* ' ,F9.3,2D13.4) 

PRINT 23, DPSIO, DTMAX, DETAO, ETALIM, DRO , RMIN, 

1 EPST , EPSR, 

2 NPHI , NTAUO ,NETA0 

23 F0RMAT(/1X, 'INTEGRATION DATA’/ 

1 IX, ’DPSIO, DTMAX, DETAO, ETALIM* ', 4F9 . 4/ 

2 IX, ’DR0,RMIN,=’,2D13.4/ 

3 IX, ’EPST,EPSR=’,2D12.3/ 

4 IX, 'NPHI, NTAUO, NETAO*’, 316) 

PRINT 24, ENA , UA , PS I A* DEG, PHIA*DEG 

24 F0RMAT(/1X, 'ABBREVIATED AIR DATA’/ 

1 IX, 'ENA,UA=',2D13.4/ 

2 IX, ’PSIA,PHIA=',2F9.1) 

GO TO (251,252), I FAN 

251 CONTINUE 

PRINT 2510, I FAN 

2510 FORMATC/1X, ’RING-FAN APPROXIMATED AS PLANAR. IFAN=',I4) 

GO TO 250 

252 CONTINUE 

PRINT 2520, IFAN 

2520 FORMATC/1X, 'RING-FAN APPROXIMATED BY MATCHED APPROXIMATION.’, 

1 4X, 'IFAN*', 14) 

250 CONTINUE 
PRINT 29 

29 FORMATC///1X, 'END DATA'///) 

I F( I FIRST . EQ . 0 . AND . I FAN . EQ . 2) 

1CALL HMSET 

IF( I FAN . EQ . 2) IFIRST = IFIRST+1 

RETURN 

END 

C$OPTIONS LIST 

SUBROUTINE SOF(ISTOP) 

IMPLICIT REAL*3(A-H,0-Z) 

CHARACTERS ISTOP(l) 

COMMON /GAMA/G, G1 ,G2,G3,G4,G5,G6,G7,G8,G9,G10,G11,G12,G13,G14,G15, 
1 G16,G17,G18,G19,G20 

COMMON /PAR/CO, ENO, EMI, D, SIGMA, TLIM, DRO, EL 0,Q0, TO, FACT, ALOGF, 

1 DPSIO, DTMAX, DETAO, ETALIM, XSI,XSF 

COMMON /NPAR/NPHI,IPAR,NP,NR,NX,NXS,NS,NSPEC,NS1,NS2, NTAUO, NET AO, 

1 NAMB, NCASE, I CASE, IFAN 

COMMON /GEOM/APF, PAI,PAI2,W,SW,CW, BETA , SBETA , CBETA , PSI 1 , SPSI 1 , 

1 CPSI1,PSIF,SPSIF,CPSIF,TPSIF,AK,SK,CK,A0,RF,XF,YF,ZF, 

2 PHISOF,PHIF,SPHIF,CPHIF,DYMIN,RMIN,XS,DIST,XO,YO,ZO, 



AMB0361 

AMB0362 

AMB0363 

AMB0364 

AMB0365 

AMB0366 

AMB0367 

AMB0368 

AMB0369 

AMB0370 

AMB0371 

AMB0372 

AMB0373 

AMB0374 

AMB0375 

AMB0376 

AMB0377 

AMB0378 

AMB0379 

AMB0380 

AMB0381 

AMB0382 

AMB0383 

AMB0384 

AMB0385 

AMB0386 

AMB0387 

AMB0388 

AMB0389 

AMB0390 

AMB0391 

AMB0392 

AMB0393 

AMB0394 

AMB0395 

AMB0396 

AMB0397 

AMB0398 

AMB0399 

AMB0400 

AMB0401 

AMB0402 

AMB0403 

AMB0404 

AMB0405 

AMB0406 

AMB0407 

AMB0408 

AMB0409 

AMB0410 

AMB0411 

AMB0412 

AMB0413 

AMB0414 

AMB0415 

AMB0416 

AMB0417 

AMB0413 

AMB0419 

AMB0420 

AMB0421 

AMB0422 

AMB0423 

AMB0424 

AMB0425 

AMB0426 

AMB0427 

AMB0428 

AMB0429 

AMB0430 

AMB0431 

AMB0432 
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1 

71 



72 

73 



74 



3 DYO,DEG,PSIN,ST1,CT1,OMEGX,OMEGY,OMEGZ,XSV(21) AMB0433 

COMMON /EPSIL/EPSETA,EPST,EPSR AMB0434 

COMMON /EXT REM/ TEXT ( 5) , ETAEXTC5) , ETAKXT (5) , PHIEXT(5) , AMB0435 

1 PSIEXT(5),EMEXT(5),FEXT(5),WEXT(5), AMB0436 

2 TMAX( 5) , ETAKMXC5 ) , ETAMAXC5 ) , PSIMAX(5) , AMB0437 

3 EMMAX(5),FMAX(5), AMB0438 

4 RFMAXC 5) , PHIFMX(5) , PH I MAX ( 5 ) , WMAX( 5) AMB0439 

COMMON /SOFPR/C, DSUMF, DSUMD, T, ETA, DETA, SUM, DSUM, SUMU, DSUMU AMB0440 

COMMON /SUMS/SUMFC 5) ,SUMD(5),SUMU2(5) AMB0441 

COMMON /COUNTS/ICONTC, ICONTT, ICNTOT, ICNTMX, IQT0TC5) , ISHADC5) AMB0442 

COMMON /SPEC/ WAV, XC( 5) , WC( 5 ) , WRC(5) , XNAMEC 5) ,QFC(5),QDC(5), AMB0443 

1 QU2C(5) , FLUXCC5) , 0MEGAC5) , FLXU2CC5) , URMSCC5) AMB0444 

PRINT 1 , ISTOP AMB0445 

FORMAT (///IX, 2H**, 2X, 30A4, 2X, 2H**,/// ) AMB0446 

PRINT 71,NS,NP,NR,NX,IC0NTC, ICONTT AMB0447 

FORMAT (IX, ' NS, NP, NR, NX, ICONTC, ICONTT= ' ,616/) AMB0448 

IF(NS.GT.NSPEC) NS=1 AMB0449 

PRINT 72, RF, PHI FXDEG, PHISOF^DEG, HXDEG, BETA^DEG AMB0450 

F0RMAT(/1X, 'RF,PHIF,PHISOF,H,BETA=',D14.5,4F10.3/) AMB0451 

PRINT 73, C, T ,TLIM, ETA AMB0452 

F0RMAT(/1X, •C,T,TLIM,ETA= , ,4D14.5/) AMB0453 

PRINT 74, DSUM, SUM, DSUMF, SUMF( NS) , SUMD(NS) , QDC( NS ) , QFC( NS) , AMB0454 

1 FLUXC( NS) , OMEGA ( NS) AMB0455 

FORMAT ( IX, 'DSUM, SUM, DSUMF, SUMF(NS) , SUMD( NS ) = ' , 5D14 . 5/ AMB0456 

1 IX, 'QDC(NS),QFC(NS),FLUXC(NS),0MEGA(NS)=',4D14.5/) AMB0457 

XX=-1 . DO AMB0458 

YY=DSQRT(XX)+1 . DO AMB0459 

STOP AMB0460 

END AMB0461 

SUBROUTINE FLUX AMB0462 

IMPLICIT REAL*8(A-H,0-Z) AMB0463 

COMMON /GAMA/G, G1,G2,G3,G4,G5,G6,G7,G8,G9,G10,G11,G12,G13,G14,G15,AMB0464 
1 G16,G17,G18,G19,G20 AMB0465 

COMMON /PAR/CO, ENO, EMI, D, SIGMA, TLIM, DRO, EL 0,Q0, TO, FACT, ALOGF, AMB0466 

1 DPSIO , DTMAX, DETAO , ETALIM, XSI , XSF AMB0467 

COMMON /NPAR/NPHI,IPAR,NP,NR,NX,NXS,NS,NSPEC,NS1,NS2,NTAU0,NETA0, AMB0468 
1 NAMB, NCASE, I CASE, I FAN AMB0469 

COMMON /GE0M/APF,PAI,PAI2,W,SW,CW,BETA,SBETA,CBETA,PSI1,SPSI1, AMB0470 

1 CPSI1,PSIF,SPSIF,CPSIF,TPSIF,AK,SK,CK,A0,RF,XF,YF,ZF,AMB0471 

2 PHISOF,PHIF,SPHIF,CPHIF, DYMIN , RMIN,XS,DIST,XO,YO,ZO, AMB0472 

3 DY0,DEG,PSIN,ST1,CT1, OMEGX, OMEGY , OMEGZ, XSV( 21 ) AMB0473 

COMMON /EPSIL/EPSETA,EPST,EPSR AMB0474 

COMMON / EXTREM/TEXT ( 5 ) , ETA EXT (5), ETAKXT (5), PHI EXT ( 5 ) , AMB0475 

1 PSI EXT ( 5 ) , EMEXT ( 5 ) , FEXT ( 5 ) , WEXT ( 5 ) , AMB0476 

2 TMAX(5),ETAKMX(5), ETAMAX( 5 ) , PSI MAX ( 5 ) , AMB0477 

3 EMMAX(5 ) , FMAX( 5) , AMB0478 

4 RFMAX( 5 ) , PHI FMX( 5),PHIMAX(5), HMAX( 5 ) AMB0479 

COMMON /SOFPR/C, DSUMF, DSUMD, T , ETA , DETA , SUM, DSUM, SUMU, DSUMU AMB0480 

COMMON /COUNTS/ ICONTC, ICONTT, ICNTOT, ICNTMX, IQTOT( 5), I SHAD( 5) AMB0481 

COMMON /SPEC/WAV, XC(5),WC( 5 ),WRC(5),XNAME(5),QFC(5),QDC( 5), AMB0482 

1 QU2C(5),FLUXC(5), OMEGA( 5) ,FLXU2C(5), URMSC( 5) AMB0483 

COMMON /SUMS/SUMF(5),SUMD(5),SUMU2(5) AMB0484 

ELO = SIGMA5(RF5<ENO AMB0435 

IF(ZO.NE.O.) AMB0486 

1CALL SOFCTHE SCHEME HERE IS NOT WRITTEN FOR ZO.NE.O.’) AMB0437 

YY0=(Y0-A0)/X0 AMB0488 

PCHECK=DATAN( YYO ) AMB0489 

I F( PCHECK . GT . PSI F-l . D-4 . OR . PCHECK . LT . -1 . D-4 ) AMB049 0 

1CALL SOF( ' FLUX RECEIVING POINT WITHIN FAN OR WITHIN SPACECRAFT') AMB0491 
SPHIF=DSIN(PHIF) AMB0492 

CPHIF=DCOS(PHIF) AMB0493 

XF=RFXCPSIF AMB0494 

YF=APFXCPHIF AMB0495 

ZF=APF^SPHIF AMB0496 

TBETA=ZF/(YF-YO) AMB0497 

BET A = DATAN( TBETA ) AMB0498 

IF(DABS(BETA) .GT.PAI2) CALL SOF( * BETA . GT . PAI/2 • ) AMB0499 

SBETA=DSIN(BETA) AMB0500 

CBETA=DCOS( BETA) AMB0501 

DIST=DSQRT( ( XF-XO )**2+( YF-YO )**2+(ZF-ZO )**2) AMB0502 

CW=(XF-XO)/DIST AMB0503 

SW=DSQRT( 1 • DO-CWXX2) AMB0504 
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H=PAI2-DATAN(CW/SW) 

OMEGX=CW 
OMEGY=SW*CBETA 
OMEGZ=SW*SBETA 
CALL LIMIT 

DO 20 NS=NS1,NS2 
SUMF(NS)=0. 

SUMU2 C NS ) =0 . 

SUMD(NS)=0 . 

FEXT (NS)=0 . 

CALL SUMT 

SUMF( NS ) =SUM 

SUMU2(NS)=SUMU 

QFC(NS)=SUMF(NS)/FACT 

QU2C(NS)=SUMU2(NS)/FACT 

FEXT(NS) =FEXT (NS )/FACT 

CALL FAN (TEXT ( NS) , PS I EXT (NS) , PHI EXT (NS) ) 

IF( PS I EXT (NS) .LT.PSIF-1 .D-10) CALL SOF( ' PS I EXT (NS).LT.PSIF') 

IF( PSIEXT (NS) .GT.PSI1) PSIEXT(NS) =PSI1 
PSIO=PSIEXT(NS) 

T =TEXT ( NS) 

CALL MATCH ( T , PSI 0 , EM, TETA) 

EMEXT(NS)=EM 
NEXT ( NS )=W 

IQTOT(NS)=IQTOT (NS)+ICONTT 
20 CONTINUE 
RETURN 
END 

SUBROUTINE LIMIT 
IMPLICIT REAL*8(A-H,0-Z) 

COMMON /GAMA/G,G1,G2,G3,G4,G5,G6,G7,G8,G9,G10,G11,G12,G13,G14,G15, 
1 G16,G17,G18,G19,G20 

COMMON /PAR/CO, ENO, EMI, D, SIGMA, TLIM,DRO, EL 0,Q0, TO, FACT, ALOGF, 

1 DPSIO,DTMAX,DETAO,ETALIM,XSI,XSF 

COMMON /NPAR/NPHI, IPAR, NP, NR, NX, NXS, NS, NSPEC, NS1 , NS2, NTAUO, NETAO , 

1 NAMB , NCASE, I CASE, I FAN 

COMMON /GE0M/APF,PAI,PAI2,W,SW,CW, BETA, SBETA,CBETA, PSI 1,SPSI1, 

CPSI1,PSIF,SPSIF,CPSIF,TPSIF,AK,SK,CK,A0,RF,XF,YF,ZF, 
PHISOF,PHIF,SPHIF,CPHIF, DYMIN, RMIN, XS , DIST, XO, YO , ZO, 
DYO,DEG,PSIN,ST1,CT1,OMEGX,OMEGY,OMEGZ,XSV(21) 

COMMON /EPSIL/EPSETA, EPST, EPSR 

COMMON /EXTREM/TEXT(5),ETAEXT(5),ETAKXT(5),PHIEXT(5), 

PSIEXT (5),EMEXT(5), FEXT ( 5) , WEXT ( 5 ) , 
TMAX(5),ETAKMX(5),ETAMAX(5),PSIMAX(5), 
EMMAX(5),FMAX(5), 

RFMAX( 5 ) , PHI FMX( 5 ) , PHIMAX( 5 ) , WMAX( 5) 

COMMON /SPEC/WAV, XC( 5) , WC(5) , WRC( 5) , XNAME( 5) , QFC(5) , QDC( 5) , 

1 QU2C(5),FLUXC(5),0MEGA(5),FLXU2C(5),URMSC(5) 

AAA=(CW/CPSI1)**2-1 . DO 
IF( AAA .LT.l.D-10) GO TO 1 
TPSI1=SPSI1/CPSI1 
AP1=A0+XF*TPSI1 

BBB = 2 . D0*( AP1XCWXTPSI1-SW*APFX(CBETAXCPHIF+SBETAXSPHIF) ) 

CCC=AP1X*2-APF**2 

DDD = BBB**2-4 . D0*AAA*CCC 

TL IM=(-BBB+DSQRT (DDD))/(2. D0*AAA) 

TLIM=TLIM/RF 

RETURN 

1 CONTINUE 

TL IM=1 . D 55 

RETURN 

END 

SUBROUTINE SUMT 
IMPLICIT REAL*8(A-H,0-Z) 

COMMON /GAMA/G,G1,G2,G3,G4,G5,G6,G7,G8,G9,G10,G11,G12,G13,G14,G15, 
1 G16,G17,G18,G19,G20 

COMMON /PAR/CO, ENO, EMI, D, SIGMA, TLIM,DR0, EL 0,Q0, TO, FACT, ALOGF, 

1 DPSIO, DTMAX, DETAO , ETAL IM, XSI , XSF 

COMMON /NPAR/NPHI, IPAR,NP, NR, NX, NXS, NS, NSPEC, NS1 , NS2, NTAUO, NETAO, 

1 NAMB, NCASE, I CASE, I FAN 

COMMON /GE0M/APF,PAI,PAI2,W,SW,CW,BETA,SBETA,CBETA,PSI1,SPSI1, 



AMB0505 
AMB0506 
AMB0507 
AMB0508 
AMB0509 
AMB0510 
AMB0511 
AMB0512 
AMB0513 
AMB0514 
AMB0515 
AMB0516 
AMB0517 
AMB0518 
AMB0519 
AMB0520 
AMB0521 
AMB0522 
AMB0523 
AMB0524 
AMB0525 
AMB 0526 
AMB0527 
AMB0528 
AMB0529 
AMB0530 
AMB 0531 
AMB0532 
AMB0533 
AMB0534 
AMB0535 
AMB 0536 
AMB0537 
AMB 0538 
AMB0539 
AMB0540 
AMB0541 
AMB 0542 
AMB0543 
AMB0544 
AMB 0 545 
AMB0546 
AMB0547 
AMB0548 
AMB0549 
AMB0550 
AMB0551 
AMB0552 
AMB0553 
AMB0554 
AMB0555 
AMB0556 
AMB 0557 
AMB0558 
AMB0559 
AMB0560 
AMB0561 
AMB 056 2 
AMB0563 
AMB0564 
AMB0565 
AMB0566 
AMB0567 
AMB0568 
AMB0569 
AMB0570 
AMB0571 
AMB0572 
AMB0573 
AMB0574 
AMB0575 
AMB0576 



41 



\MB 



SCRIPT A1 



1 CPSI1,PSIF,SPSIF,CPSIF,TPSIF,AK,SK,CK,A0,RF,XF,YF,ZF,AMB0577 

2 PHISOF,PHIF,SPHIF,CPHIF,DYMIN,RMIN,XS,DIST,XO,YO,ZO, AMB0578 

3 DYO,DEG,PSIN,ST1,CT1,OMEGX,OMEGY,OMEGZ,XSVC21) AMB0579 

COMMON / EPS I L/ EPS ETA, EPST , EPSR AMB0580 

COMMON /EXT REM/ TEXT C5 ) , ETA EXT C 5) , ETAKXTC5) , PHI EXT C 5) , AMB0581 

1 PSIEXT C 5) , EMEXT C 5) , FEXT C5 ) , WEXT C 5) , AMB0582 

2 TMAXC 5) , ETAKMXC 5) , ETAMAXC5) ,PSIMAXC5), AMB0583 

3 EMMAXC5),FMAXC5), AMB0584 

4 RFMAXC 5) , PHIFMXC 5) , PHIMAX(5) , WMAXC 5) AMB0585 

COMMON /SOFPR/CC, DSUMF, DSUMD, T , ETA , DETA, SUM, DSUM, SUMU, DSUMU AMB0586 

COMMON /COUNTS/ ICO NTC, ICO NTT, ICNTOT , ICNTMX, IQT0T(5),ISHAD(5) AMB0587 

COMMON /SPEC/WAV, XC( 5),WC(5),WRC(5),XNAME(5),QFC(5),QDC(5), AMB0588 

1 QU2CC5) , FLUXCC5) , OMEGA (5) ,FLXU2C(5), URMSCC5) AMB0589 

COMMON /SUMS/SUMFC5 ) , SUMD( 5) , SUMU2( 5) AMB0590 

: INTEGRATION OF FLUX ARRIVING ALONG A SINGLE RAY AMB0591 

DT=DPSI0 AMB0592 

PSIN=PSIF AMB0593 

ETA1=0 . AMB0594 

ETA3=0 . AMB0595 

FETA4=0 . AMB0596 

FETAU4=0. AMB0597 

T = 0 . AMB0598 

SUM=0 . AMB0599 

SUMU=0. AMB0600 

ICONTT=0 AMB0601 

1 ICONTT=ICONTT+l AMB0602 

PSIL=PSIN AMB0603 

DT2=DT/2 . DO AMB0604 

DT6=DT/6 . DO AMB0605 

T1=T+DT2 AMB0606 

T2=T+DT AMB0607 

FETA1=FETA4 AMB0608 

FETAU1=FETAU4 AMB0609 

CALL PATHKCT1 , ETAK1 ) AMB0610 

CALL FETA( T1 , ETA1 , ETAK1 , GT2, FETA2, FETAU2) AMB0611 

FETA3=FETA2 AMB0612 

FETAU3=FETAU2 AMB0613 

CALL PATHKCT2, ETAK3) AMB0614 

CALL FETA(T2, ETA3, ETAK3, GT4, FETA4, FETAU4) AMB0615 

DETA=DT*GT4 AMB0616 

DSUM=DT65(( FETA 1+2 . D0*( FETA2 + FETA3 ) + FETA4) AMB0617 

DSUMU=DT6X(FETAU1+2.D0X(FETAU2+FETAU3)+FETAU4) AMB0618 

T=T+DT AMB0619 

ETA=ETA3 AMB0620 

ETAK=ETAK3 AMB0621 

SUM=SUM+DSUM AMB0622 

SUMU =SUMU+ DSUMU AMB0623 

IF(FEXT(NS) .GT.FETA4) GO TO 10 AMB0624 

FEXT(NS)=FETA4 AMB0625 

TEXT C NS ) =T AMB0626 

ETA EXT (NS)=ETA AMB0627 

ETAKXT (NS)=ETAK AMB0628 

10 CONTINUE AMB0629 

: STEP CONTROL (DT) AMB0630 

CALL FAN( T , PSI , PHI ) AMB0631 

IFCPSI.LT. PSIF-l.D-10) CALL SOFC 'PSI.LT.PSIF') AMB0632 

IFCPSI.GT.PSIl) PSI=PSI1 AMB0633 

PSIN=PSI AMB0634 

DPSI=PSIN-PSIL AMB0635 

DTP=DT^(DPSIO/(DPSI+1.D-10)) AMB0636 

DTE=DTx( DETAO/ ( DETA+1 . D-10 ) ) AMB0637 

DT1=1 .2D0*DT AMB0638 

DT=DMIN1(DTP, DTE, DTI, DTMAX) AMB0639 

IFCDT.LE.O.) CALL SOFC ' COMPUTED DT NEGATIVE') AMB0640 

15 CONTINUE AMB0641 

IFCIPAR.LT. 1) AMB0642 

1PRINT 111,NR,NP,T,PSI*DEG,PHI*DEG,ETA, ETAK, SUM, DSUM/ CSUM+1 .D-20) AMB0643 

111 FORMAT C1X,'NR,NP,T,PSI,PHI = , ,2I3,3D12.3/ AMB0644 

1 IX, ’ETA,ETAK,SUM,ERRR=',4D12.3) AMB0645 

IFCICONTT.GT.ICNTMX) AMB0646 

1CALL SOFC ' ICONTT TOO LARGE') AMB0647 

IFC ICONTT . LE . 2) GO TO 1 AMB0648 
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I FC ETA+ETAK . GT . ETALIM) GO TO 100 AMB0649 

IF(T . GT . 50 . DO .OR. TXRF.GT.AO) GO TO 100 AMB0650 

IF(SUM. EQ .0 . ) GO TO 1 AMB0651 

ERR=CDSUM/SUM)/DT AMB0652 

IF(ERR.GT.EPST) GO TO 1 AMB0653 

100 CONTINUE AMB0654 

SUM=SUM*ELO AMB0655 

SUMU=SUMU*ELO AMB0656 

RETURN AMB0657 

END AMB0658 

SUBROUTINE FETACT, ETAIK, ETAK, GT, FET, FETU2) AMB0659 

IMPLICIT REAL*3CA-H,0-Z) AMB0660 

REALXS MU1,MU2 AMB0661 

COMMON /GAMA/G,G1,G2,G3,G4,G5,G6,G7,G8,G9,G10,G11,G12,G13,G14,G15, AMB0662 

1 G16,G17,G18,G19,G20 AMB0663 

COMMON /PAR/CO, ENO, EMI, D, SIGMA, TLIM,DRO, EL 0,Q0, TO, FACT, ALOGF, AMB0664 

1 DPSIO , DTMAX, DETAO, ETALIM, XSI , XSF AMB0665 

COMMON /NPAR/NPHI,IPAR,NP,NR,NX,NXS,NS,NSPEC,NS1,NS2,NTAU0,NETA0, AMB0666 

1 NAMB, NCASE, ICASE, I FAN AMB0667 

COMMON /GEOM/APF, PAI , PAI2,W,SW,CW,BETA,SBETA,CBETA,PSI1,SPSI1, AMB0668 

1 CPSI1,PSIF,SPSIF,CPSIF,TPSIF,AK,SK,CK,A0,RF,XF,YF,ZF,AMB0669 

2 PHISOF,PHIF,SPHIF, CPHIF, DYMIN, RMIN,XS,DIST,X0,Y0,Z0, AMB0670 

3 DY0,DEG,PSIN,ST1,CT1, OMEGX, OMEGY , OMEGZ, XSV ( 21 ) AMB0671 

COMMON /EPSIL/EPSETA, EPST, EPSR AMB0672 

COMMON /EXT REM/ TEXT ( 5) , ETA EXT ( 5) , ET AKXT C 5) , PH I EXT ( 5) , AMB0673 

1 PSIEXT ( 5) , EM EXT ( 5) , FEXT ( 5) , WEXT ( 5) , AMB0674 

2 TMAX( 5 ) , ETAKMXC 5) , ETAMAXC 5) , PSIMAXC 5) , AMB0675 

3 EMMAXC 5) , FMAX(5) , AMB0676 

4 RFMAXC 5) , PHIFMXC 5) , PHIMAX(5) , WMAXC 5) AMB0677 

COMMON /SPEC/WAV ,XC( 5 ),WCC5),WRCC5),XNAMEC5),QFCC5),QDCC5), AMB0678 

1 QU2CC 5) , FLUXC( 5) , OMEGA ( 5) , FLXU2CC 5) , URMSC( 5 ) AMB0679 

COMMON /AMBIEN/ENA, UA,PSIA,PHIA,HA(3),WA(3), AMB0680 

1 UAX, UAY, UAZ, AA, BA, CA, RA, XA, YA, ZA, SHADOW AMB0681 

LOGICAL SHADOW AMB0682 

COMMON /NAGESH/PIK, UIK,UIKX,UIKY,UIKZ AMB0683 

ETAIK=0. AMB0684 

IF( SHADOW) GO TO 1 AMB0685 

K= 1 AMB0636 

I=NS AMB0637 

CALL FAN( T , PSI , PHI ) AMB0638 

IFCPSI.LT. PSIF-l.D-10) CALL SOFC 'PSI.LT.PS IF') AMB0639 

I FCPSI . GT . PSI1 ) PSI = PSI1 AMB0690 

PSIO=PSI AMB0691 

CALL MATCHCT , PSIO , EM, TETA) AMB0692 

SPSI=DSIN( PSI ) AMB0693 

CPSI = DCOS( PSI ) AMB0694 

SPHI = DSIN( PHI ) AMB0695 

CPHI = DCOS( PHI ) AMB0696 

ST=DSIN(TETA) AMB0697 

CT=DCOS( TETA) AMB0698 

G0REM = 1 .D0 + G1XEMX5«2 AMB0699 

TERMN=G0REMXXG6 AMB0700 

U=EM*CO/DSQRT ( GOREM) AMB0701 

UX=UXCT AMB0702 

UY=UXSTXCPHI AMB0703 

UZ=UXSTXSPHI AMB0704 

: COLLISION AMB0705 

MU1=WC(I)/(WC(I)+WA(K)) AMB0706 

MU2=1 . D0-MU1 AMB0707 

UMX=MU1XUX+MU2XUAX AMB0708 

UMY=MU1XUY+MU2XUAY AMB0709 

UMZ=MU1XUZ+MU2XUAZ AMB0710 

D0TUM=0MEGXXUMX+0MEGYXUMY+0MEGZXUMZ AMB0711 

URX=UX-UAX AMB0712 

URY=UY-UAY AMB0713 

URZ=UZ-UAZ AMB0714 

UR = DSQRTCURXX5«2+URYX5(2+URZX5«2) AMB0715 

DET = D0TUMXX2+CMU2XUR)X5(2-CUMXXX2+UMYX5(2+UMZX5(2) AMB0716 

IFCDET.LT. 0.) GO TO 1 AMB0717 

DET1 = DSQRT C DET ) AMB0718 

UIK1=-D0TUM+DET1 AMB0719 

UIK2=-D0TUM-DET1 AMB0720 
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IFCUIK2.GT.0.) CALL SOFC'DOUBLE COLLISION OPTION NOT PROGRAMMED AMB0721 
1 YET') AMB0722 

UIK=UIK1 AMB0723 

IFCUIK.LE.O. ) GO TO 1 AMB0724 

UIKX=-OMEGX*UIK AMB0725 

UIKY=-OMEGY*UIK AMB0726 

UIKZ=-OMEGZ*UIK AMB0727 

CDEL=CD0TUM+UIK)/CMU2*UR) AMB0728 

IFCCDEL . LE. 0 . ) CALL SOFC*CDEL NEGATIVE NOT PROGRAMMED YET*) AMB0729 

IFCCDEL-1 .D-10.GT.1 .DO) AMB0730 

1CALL SOFC'CDEL (COS( DELTA) ) CANNOT BE GT.l. 1 ) AMB0731 

PIK=CUIK/CMU2*UR))XX2/C4.D0XPAIXCDEL) AMB0732 

IF CPIK.LT.O.) CALL SOF( • PIK . LT . 0 1 ) AMB0733 

FET=CUR/UA)XPIK/TERMN AMB0734 

UREL = DSQRTC CUX-UIKX)XX2+CUY-UIKY)XX2+CUZ-UIKZ)**2) AMB0735 

GT=ELOX(UREL/UIK)/TERMN AMB0736 

CALL PATHIKC T , ETAIK) AMB0737 

POWER=ETAIK+ETAK-ALOGF AMB0738 

EFACT=0 . AMB0739 

I F( POWER . LT . 6 0 . DO) EFACT = DEXP( -POWER) AMB0740 

FET=FETXEFACT AMB0741 

FETU2=FETXUIKXX2 AMB0742 

IFCEM.LT. 0.) CALL SOFC * EM . LT . 0 * ) AMB0743 

FETU2=FET*EM AMB0744 

RETURN AMB0745 

CONTINUE AMB0746 

FET=0 . AMB0747 

FETU2=0. AMB0748 

GT=0. AMB0749 

RETURN AMB0750 

END AMB0751 

SUBROUTINE PATHIKC TC, ETAIK) AMB0752 

IMPLICIT REALX8(A-H,0-Z) AMB0753 

REAL*8 MU1,MU2 AMB0754 

COMMON /GAMA/G,G1,G2,G3,G4,G5,G6,G7,G8,G9,G10,G11,G12,G13,G14,G15,AMB0755 
1 G16,G17,G18,G19,G20 AMB0756 

COMMON /PAR/CO, ENO, EMI, D, SIGMA, TLIM,DR0, EL 0,Q0, TO, FACT, ALOGF, AMB0757 

1 DPSIO,DTMAX,DETAO,ETALIM,XSI,XSF AMB0758 

COMMON /NPAR/NPHI,IPAR,NP,NR,NX,NXS,NS,NSPEC,NS1,NS2,NTAU0, NETAO, AMB0759 
1 NAMB , NCASE, ICASE, I FAN AMB0760 

COMMON /GE0M/APF,PAI,PAI2, W,SW,CW,BETA,SBETA,CBETA,PSI1,SPSI1, AMB0761 
CPSI1,PSIF,SPSIF,CPSIF,TPSIF,AK,SK,CK,A0,RF,XF,YF,ZF,AMB0762 
PHISOF,PHIF,SPHIF,CPHIF, DYMIN, RMIN, XS , DIST, XO , YO , ZO , AMB0763 
DYO , DEG, PSIN, ST1 ,CT1 , OMEGX, OMEGY, OMEGZ, XSVC 21 ) AMB0764 

COMMON /EPSIL/EPSETA, EPST, EPSR AMB0765 

COMMON /EXT REM/ TEXT C5),ETAEXTC5), ETAKXT C 5 ) , PH I EXT C 5) , AMB0766 

P SI EXT C 5) , EM EXT C 5 ) , FEXT C5),WEXTC5), AMB0767 

TMAXC 5 ) , ETAKMXC 5),ETAMAX(5),PSIMAXC5), AMB0768 

EMMAXC 5 ) , FMAXC 5 ) , AMB0769 

RFMAXC 5 ) , PHIFMXC 5) , PHIMAXC5), WMAXC 5) AMB077 0 

COMMON /SPEC/ WAV , XCC 5 ) , WCC 5),WRC(5), XNAMEC 5),QFC(5),QDCC5), AMB0771 

1 QU2CC 5 ) , FLUXCC 5 ) , OMEGA C 5 ),FLXU2CC5), URMSCC 5) AMB0772 

COMMON /AMB IEN/ENA,UA,PSIA,PHIA,HAC3),WAC3), AMB0773 

1 UAX, UAY , UAZ, AA, BA, CA, RA, XA, YA, ZA, SHADOW AMB0774 

LOGICAL SHADOW AMB0775 

N ETA = NETAO AMB0776 

DT=TC/DBLE(NETA) AMB0777 

DT2=DT/2 . DO AMB0778 

DT6=DT/6 . DO AMB0779 

GT4=0 . AMB0780 

T=0. AMB0781 

ETA=0 . AMB0782 

IT=0 AMB0783 

IT=IT+1 AMB0784 

T1=T+DT2 AMB0785 

T2=T+DT AMB0786 

GT1 =GT 4 AMB0787 

CALL FTC T1 , GT2) AMB0788 

GT3=GT2 AMB0789 

CALL FTC T2, GT4) AMB0790 

DETA = DT6X( GT1+2 . DO*C GT2+GT3)+GT4) AMB0791 

T=T+DT AMB0792 



44 



AMB 



SCRIPT A1 



ETA=ETA+DETA AMB0793 

IFCIT.LT. NETA) GO TO 1 AMB0794 

I™ IK* ETA AMB0795 

RETURN AMB0796 

END AMB0797 

SUBROUTINE FTCT,GT) AMB0798 

IMPLICIT REALXSC A-H, O-Z) AMB0799 

REAL*8 MU1,MU2 AMB0800 

COMMON /GAMA/G,G1,G2,G3,G4,G5,G6,G7,G8,G9,G10,G11,G12,G13,G14,G15,AMB0801 
1 G16,G17,G18,G19,G20 AMB0802 

COMMON /PAR/CO, ENO, EMI, D, SIGMA, TLIM, DRO , ELO , QO , TO, FACT, ALOGF, AMB0803 

1 DPSIO,DTMAX,DETAO,ETALIM,XSI,XSF AMB0804 

COMMON /NPAR/NPHI,IPAR,NP,NR,NX,NXS,NS,NSPEC,NS1,NS2,NTAU0,NETA0, AMB0805 
1 NAMB , NCASE, I CASE, I FAN AMB0806 

COMMON /GE0M/APF,PAI,PAI2,W,SW,CW,BETA,SBETA,CBETA,PSI1,SPSI1, AMB0807 

1 CPSI1,PSIF,SPSIF,CPSIF,TPSIF,AK,SK,CK,A0,RF,XF,YF,ZF,AMB0808 

2 PHISOF,PHIF,SPHIF,CPHIF,DYMIN,RMIN,XS,DIST,XO,YO,ZO, AMB0809 

3 DYO,DEG,PSIN,ST1,CT1,OMEGX,OMEGY,OMEGZ,XSVC21) AMB0810 

COMMON /EPSIL/EPSETA,EPST,EPSR AMB0811 

COMMON /EXTREM/TEXTC5),ETAEXTC5),ETAKXTC5),PHIEXTC5), AMB0812 

1 PS I EXT C 5) , EM EXT C 5) , FEXT C 5) , WEXT C5 ) , AMB0813 

2 TMAXC5),ETAKMXC5),ETAMAXC5),PSIMAXC5), AMB0814 

3 EMMAXC 5 ) , FMAXC 5) , AMB0815 

4 RFMAXC5) , PHI FMXC 5),PHIMAXC5) , WMAXC 5) AMB0816 

COMMON /SPEC/WAV,XCC5),WCC5),WRCC5),XNAMEC5),QFCC5),QDCC5), AMB0817 

1 QU2CC 5 ) , FLUXCC5) , OMEGAC 5) ,FLXU2CC5) , URMSCC 5) AMB0818 

COMMON /AMBIEN/ENA,UA,PSIA,PHIA,HAC3),WA(3), AMB0819 

1 UAX,UAY,UAZ,AA, BA, CA , RA , XA, YA, ZA , SHADOW AMB0820 

LOGICAL SHADOW AMB0821 

COMMON /NAGESH/PIK,UIK,UIKX,UIKY,UIKZ AMB0822 

K=1 AMB0823 

I=NS AMB0824 

CALL FANCT,PSI,PHI) AMB0825 

IFCPSI .LT.PSIF-1 .D-10) CALL SOFC ' PSI . LT . PSIF* ) AMB0826 

IFCPSI .GT.PSI1) PSI=PSI1 AMB0827 

PSI0=PSI AMB0828 

CALL MAT CHCT,PSIO,EM,TETA) AMB0829 

SPSI=DSIN(PSI) AMB0830 

CPSI = DCOS(PSI ) AMB0831 

SPHI =DSIN( PHI ) AMB0832 

CPHI = DCOS ( PHI ) AMB0833 

ST=DSINCTETA) AMB0834 

CT=DCOSCTETA) AMB0835 

GOREM=l .D0+G1*EM**2 AMB0836 

TERMN=GOREM**G6 AMB0837 

U=EM*CO/DSQRTCGOREM) AMB0833 

UX=U*CT AMB0839 

UY=U*ST*CPHI AMB0840 

UZ=U*ST*SPHI AMB0841 

UREL=DSQRTCCUX-UIKX)**2+CUY-UIKY)X*2+CUZ-UIKZ)X*2) AMB0842 

GT=ELO*CUREL/UIK)/TERMN AMB0843 

RETURN AMB0844 

END AMB0845 

SUBROUTINE PATHKCT , ETAK) AMB0846 

IMPLICIT REAL*8CA-H,0-Z) AMB0847 

COMMON /GAMA/G,G1,G2,G3,G4,G5,G6,G7,G8,G9,G10,G11,G12,G13,G14,G15,AMB0848 
1 G16,G17,G18,G19,G20 AMB0849 

COMMON /PAR/CO, ENO, EMI, D, SIGMA, TLIM, DRO, EL 0,Q0, TO, FACT, ALOGF, AMB0850 

1 DPSIO,DTMAX, DETAO, ETALIM,XSI,XSF AMB0851 

COMMON /NPAR/NPHI, IPAR , NP , NR , NX, NXS, NS , NSPEC, NS1 , NS2, NTAUO , NETAO , AMB0852 
1 NAMB, NCASE, I CASE, IFAN AMB0853 

COMMON /GEOM/APF ,PAI,PAI2,W,SW,CW,BETA,SBETA,CBETA,PSI1,SPSI1, AMB0854 

1 CPSI1,PSIF,SPSIF,CPSIF,TPSIF,AK,SK,CK,A0,RF,XF,YF,ZF,AMB0855 

2 PHISOF,PHIF,SPHIF,CPHIF,DYMIN,RMIN,XS,DIST,XO,YO,ZO, AMB0856 

3 DY0,DEG,PSIN,ST1,CT1, OMEGX, OMEGY , OMEGZ, XSV ( 21 ) AMB0857 

COMMON /EPSIL/EPSETA, EPST, EPSR AMB0858 

COMMON /EXT REM/ TEXT C5),ETAEXTC5),ET AKXT ( 5 ) , PHI EXT ( 5) , AMB0859 

1 PS I EXT (5) , EM EXT C5),FEXTC5), WEXT C5 ) , AMB0860 

2 TMAXC 5 ) , ETAKMXC 5) ,ETAMAX(5),PSIMAX(5), AMB0861 

3 EMMAXC 5 ) , FMAXC 5) > AMB0862 

4 R FMAXC 5 ) , PH I FMXC 5) ,PHIMAXC5), WMAXC 5) AMB0863 

COMMON /COUNTS/ICONTC, ICONTT, ICNTOT , ICNTMX, IQTOTC 5) , ISHADC 5) AMB0864 
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COMMON /AMBIEN/ENA,UA,PSIA,PHIA,HAC3),WAC3), AMB0865 

1 UAX,UAY,UAZ,AA, BA, CA, RA, XA, YA, ZA, SHADOW AMB0866 

LOGICAL SHADOW AMB0867 

ETAK=0 . AMB0868 

DETERMINE POINT OF ENTRY OF AMBIENT TRAJECTORY TO FAN AMB0869 

TRF=T*RF AMB0870 

XC=XF+TRF*OMEGX AMB0871 

YC=YF+TRF*OMEGY AMB0872 

ZC=ZF+TRF*OMEGZ AMB0873 

CHECK SHADOH AMB0874 

SHADOH=. FALSE. AMB0875 

EVER=BAX*2+CAXX2 AMB0876 

DETS=EVER*A0*X2-CBAXZC-CAXYC)X*2 AMB0877 

IFCDETS . LE. 0 . ) GO TO 2 AMB0878 

DETS1 =DSQRTC DETS) AMB0879 

TAU1=C-CBA*YC+CA*ZC)+DETS1)/EVER AMB0880 

IFCTAU1 . GT . 0 . ) SHADOW 3 . TRUE . AMB0881 

CONTINUE AMB0882 

IFCSHADOW) GO TO 10 AMB0883 

EVER1=A0+XC*TPSIF AMB0884 

EVER2=BA*x2+CAxx2-CAA*TPSIF)*x2 AMB0885 

EVER3=BAXYC+CAXZC-AAXEVER1XTPSIF AMB0886 

DET=EVER3**2-EVER2*CYC**2+ZC*x2-EVERl**2) AMB 0887 

IFCDET.LE.O.) AMB0888 

1CALL SOFC * NO INTERSECTION OF AMB. TRAJ . WITH LIMITING CONE*) AMB0889 

DET1=DSQRT CDET) AMB0890 

TAUP=(-EVER3+DET1)/EVER2 AMB0891 

TAUM=C-EVER3-DET1)/EVER2 AMB0892 

IFCTAUP .GT . 0 . .AND. TAUM.GT.O.) AMB0893 

1CALL SOFC 'TWO POSITIVE INTERSECTIONS WITH LIMITING CONE. NOT PERMITAMB0894 



DO 

DO 



10 



1 IN THIS VERSION') 

TAUF=DMAX1 CTAUP, TAUM) 

IFITAUF. LE. 0 . ) 

1CALL SOFC ' NO POSITIVE INTERSECTION WITH LIMITING CONE') 
XA=XC+TAUF*AA 
YA=YC+TAUFXBA 
ZA=ZC+TAUFXCA 

RA=DSQRT(XAXX2+(DSQRT(YAXX2+ZAXX2)-A0)XX2) 

TAUF=TAUF/RA 

NTAU=NTAUO 

DTAU=TAUF/DBL EC NTAU) 

ETAK=0 . 

TAU=0. 

DTAU2=DTAU/2, 

DTAU6 =DTAU/6 , 

GTAU4=0. 

ITAU=0 
ITAU=ITAU+1 
TAU1=TAU+DTAU2 
TAU2=TAU+DTAU 
GTAU1 =GTAU4 
CALL FTAUC TAU1 , GTAU2) 

GTAU3=GTAU2 
CALL FTAUCTAU2, GTAU4) 

DETAK=DTAU6XCGTAU1+2.D0XCGTAU2+GTAU3)+GTAU4) 
TAU=TAU+DTAU 
ETAK=ETAK+DETAK 
IFCITAU.LT. NTAU) GO TO 1 
ETAK=ETAKXC SIGMAXENOXRA) 

RETURN 
CONTINUE 

ISHADCNS)=ISHADCNS)+1 
RETURN 
END 

SUBROUTINE FTAUCTAU, GTAU) 

IMPLICIT REAL*8CA-H,0-Z) 



AMB0895 
AMB0896 
AMB0897 
AMB0898 
AMB0899 
AMB0900 
AMB0901 
AMB0902 
AMB0903 
AMB0904 
AMB0905 
AMB0906 
AMB0907 
AMB0908 
AMB0909 
AMB0910 
AMB0911 
AMB0912 
AMB0913 
AMB0914 
AMB0915 
AMB0916 
AMB0917 
AMB0918 
AMB0919 
AMB0920 
AMB0921 
AMB0922 
AMB0923 
AMB0924 
AMB0925 
AMB0926 
AMB0927 
AMB0928 
AMB0929 
AMB0930 

COMMON /GAMA/G,G1,G2,G3,G4,G5,G6,G7,G8,G9,G10,G11,G12,G13,G14,G15,AMB0931 
1 G16,G17,G18,G19,G20 AMB0932 

COMMON /PAR/CO, ENO, EMI, D, SIGMA, TLIM,DR0, EL 0,Q0, TO, FACT, ALOGF, AMB0933 

1 DPSIO , DTMAX, DETAO , ETALIM, XSI , XSF AMB0934 

COMMON /NPAR/NPHI, IPAR, NP, NR, NX, NXS, NS, NSPEC, NS1 , NS2, NTAUO, NETAO, AMB0935 
1 NAMB, NCASE, ICASE, I FAN AMB0936 
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C 

c 

c 

c 



COMMON /GE0M/APF,PAI,PAI2,W,SW,CW,BETA,SBETA,CBETA,PSI1,SPSI1, AMB0937 

1 CPSI 1, PSI F, SPSI F, CPSI F,TPSIF,AK,SK,CK, AO, RF,XF,YF,ZF, AMB0938 

2 PHISOF,PHIF,SPHIF,CPHIF,DYMIN,RMIN,XS,DIST,XO,YO,ZO, AMB0939 

3 DYO, DEG, PSIN, ST1,CT1, OMEGX, OMEGY, OMEGZ, XSVC21) AMB0940 

COMMON /EPSIL/EPSETA, EPST, EPSR AMB09A1 

COMMON /EXTREM/TEXT(5),ETAEXT(5),ETAKXT(5),PHIEXT(5), AMB0942 

1 PSIEXT(5),EMEXT(5),FEXT(5),WEXT(5), AMB0943 

2 TMAX(5),ETAKMX(5),ETAMAX(5),PSIMAX(5), AMB0944 

3 EMMAX(5),FMAX(5), AMB0945 

A RFMAX(5),PHIFMX(5),PHIMAX(5),WMAX(5) AMB0946 

COMMON /SPEC/WAV,XC(5),WC(5),WRC(5),XNAME(5),QFC(5),QDC(5), AMB0947 

1 QU2C(5),FLUXC(5),0MEGA(5),FLXU2C(5),URMSC(5) AMB09A8 

COMMON /AMBIEN/ENA,UA,PSIA,PHIA,HAC3),WA(3), AMB0949 

1 UAX,UAY,UAZ,AA, BA, CA,RA,XA,YA,ZA, SHADOW AMB0950 

LOGICAL SHADOW AMB0951 

CALL FANT(TAU,PSI,PHI) AMB0952 

IF(PSI . LT.PSIF-1 .D-10) CALL SOFC'PSI.LT.PSIF') AMB0953 

IFCPSI .GT.PSI1) PSI=PSI1 AMB0954 

PSIO=PSI AMB0955 

CALL MATCH (T ,PSIO,EM,TETA) AMB0956 

SPSI=DSIN(PSI) AMB0957 

CPSI =DCOS( PSI ) AMB0958 

SPHI =DSI N( PHI ) AMB0959 

CPHI=DCOS(PHI) AMB0960 

ST=DSIN(TETA) AMB0961 

CT=DCOS(TETA) AMB0962 

GOREM=l .D0+G1*EM**2 AMB0963 

TERMN=GOREM**G6 AMB096A 

U=EM*CO/DSQRT(GOREM) AMB0965 

UREL=DSQRT((CTXU-UAX)XX2+(STXCPHIXU-UAY)XX2+(STXSPHIXU-UAZ)XX2) AMB0966 

GTAU=UREL/(UAXTERMN) AMB0967 

RETURN AMB0968 

END AMB0969 

SUBROUTINE FAN( T , PSI , PHI) AMB0970 

IMPLICIT REAL*8(A-H,0-Z) AMB0971 

COMMON /PAR/CO, ENO, EMI, D, SIGMA, TLIM,DR0, EL 0,Q0, TO, FACT, ALOGF, AMB097 2 

1 DPSIO , DTMAX, DETAO,ETALIM,XSI,XSF AMB0973 

COMMON /GE0M/APF,PAI,PAI2,W,SW,CW, BETA, SBETA,CBETA, PSI 1,SPSI1, AMB097A 

1 CPSI1,PSIF,SPSIF,CPSIF,TPSIF,AK,SK,CK,A0,RF,XF,YF,ZF,AMB0975 

2 PHISOF,PHIF,SPHIF,CPHIF,DYMIN,RMIN,XS,DIST,XO,YO,ZO, AMB0976 

3 DYO , DEG, PSIN, ST1 ,CT1 , OMEGX, OMEGY ,0MEGZ,XSV(21) AMB0977 

COMMON /POINT/ XP,YP,XC OR, YCOR AMB0978 

RING FAN GEOMETRY. FAN CORNER IS AT ( 0 , AO*COS( PHI ),A0*SIN(PHI)) . AMB0979 

RF — RADIAL DISTANCE ON LIMITING CHARACTERISTIC OF POINT OF AMB0980 

ENTRANCE OF RAY. AMB0981 

DIRECTION COSINES OF RAY: OMEGX, OMEGY, OMEGZ AMB0982 

TRF=T*RF AMB0983 

X=XF+TRF*OMEGX AMB0984 

Y=YF+TRFXOMEGY AMB0985 

Z=ZF+TRF*OMEGZ AMB0986 

DY=DSQRT(YXY+ZXZ)-AO AMB0987 

IFCDABS(DY).LE.l.D-lOXAO) DY=1.D-10*A0 AMB0988 

IFCDY.LT.O. ) AMB0989 

1CALL SOF( ’POINT X,Y,X CANNOT BE CLOSER TO X-AXIS THAN RADIUS AO') AMB0990 
YY=X/DY AMB0991 

PS I = PA 1 2- DAT AN (YY) AMB0992 

PHI=DATAN(Z/Y) AMB0993 

XP=XCOR+X AMB0994 

YP=AO+DY AMB0995 

RETURN AMB0996 

END AMB0997 

SUBROUTINE FANT ( TAU , PSI , PHI ) AMB0998 

IMPLICIT REAL*8(A-H,0-Z) AMB0999 

COMMON /PAR/CO, ENO, EMI, D, SIGMA, TLIM, DRO , EL 0 , QO , TO , FACT , ALOGF, AMBIOOO 

1 DPSIO, DTMAX, DETAO,ETALIM,XSI,XSF AMB1001 

COMMON /GEOM/APF, PAI,PAI2,W,SW,CW, BETA , SBETA, CBETA, PSI 1 , SPSI 1 , AMB1002 

1 CPSI 1, PSI F, SPSI F, CPSI F, TPSI F,AK,SK,CK, AO, RF, XF, YF, ZF, AMB1 003 

2 PHISOF,PHIF,SPHIF,CPHIF, DYMIN, RMIN,XS,DIST,X0,Y0,Z0, AMB100A 

3 DYO , DEG, PSIN, ST1 , CT1 , OMEGX, OMEGY, OMEGZ, XSV ( 21 ) AMB1005 

COMMON /AMBI EN/ENA, UA ,PSIA,PHIA,HA(3),WA(3), AMB1006 

1 UAX, UAY, UAZ, AA , BA , CA, RA, XA, YA , ZA, SHADOW AMB1007 

COMMON /POINT/XP,YP,XCOR,YCOR AMB1008 
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C 

C 

C 

C 



C 



C 



LOGICAL SHADOW AMB1009 

RING FAN GEOMETRY. FAN CORNER IS AT ( 0 , AO*COS(PHI ) , A0*SIN( PHI ) ) . AMB1010 

RA — RADIAL DISTANCE ON LIMITING CHARACTERISTIC OF POINT OF AMB1011 

ENTRANCE OF RAY. AMB1012 

DIRECTION COSINES OF RAY: -AA,-BA,-CA AMB1013 

TRA=TAU*RA AMB1014 

X=XA-TRAXAA AMB1015 

Y=YA-TRAXBA AMB1016 

Z=ZA-TRAXCA AMB1017 

DY=DSQRT(YXY+Z*Z)-AO AMB1018 

IF(DABSCDY) . LE . 1 . D-10XA0) DY=1.D-10XA0 AMB1019 

IF(DY . LT . 0 . ) AMB1020 

1CALL SOF( 'POINT X,Y,X CANNOT BE CLOSER TO X-AXIS THAN RADIUS AO') AMB1021 
YY=X/DY AMB1022 

PSI=PAI2-DATAN(YY) AMB1023 

PHI=DATAN(Z/Y) AMB1024 

XP=XCOR+X AMB1025 

YP=A0+DY AMB1026 

RETURN AMB1027 

END AMB1028 

SUBROUTINE HMSET AMB1029 

SUBROUTINE NUMBER 20 AMB1030 

IMPLICIT REAL*8(A-H,0-Z,$) AMB1031 

REALX8 KAPAOB,MHINV,MINVO,M,MF,M1,M2,M3,NORM,MEXIT,LAMDOB AMB1032 

COMMON /GAMA/G,G1,G2,G3,G4,G5,G6,G7,G8,G9,G10,G11,G12,G13,G14,G15,AMB1033 
1 G16,G17,G18,G19,G20 AMB1034 

COMMON /PAR/CO, ENO, EMI, D, SIGMA, TLIM,DR0, EL0,Q0, TO, FACT, ALOGF, AMB1035 

1 DPSIO,DTMAX,DETAO, ETALIM,XSI,XSF AMB1036 

COMMON /GE0M/APF,PAI,PAI2,W,SW,CW,BETA,SBETA,CBETA,PSI1,SPSI1, AMB1037 

1 CPSI1,PSIF,SPSIF,CPSIF,TPSIF,AK,SK,CK,A0,RF,XF,YF,ZF,AMB1038 

2 PHISOF,PHIF,SPHIF,CPHIF,DYMIN,RMIN,XS,DIST,XO,YO,ZO, AMB1039 

3 DYO,DEG,PSIN,ST1,CT1,OMEGX,OMEGY,OMEGZ,XSV(21) AMB1040 

COMMON /GRP/DMINV,MHINV(101),HMV(101) AMB10A1 

COMMON /IGRP/KHM AMB1042 

A ROUTINE FOR THE C+ DERIVATIVE DUE TO RING SYMMETRY (GRP). AMB1043 

MEXIT = EM1 AMBlO^jA 

KHM=51 AMB1045 

IFCKHM . GT . 101 ) CALL SOF('2001') AMB1046 

MINV0=1 .DO/MEXIT AMB1047 

DMINV=MINVO/DBLE( KHM-1 ) AMB1048 

M=MEXIT AMB1049 

SUM=0 . AMB1050 

KHM1=KHM-1 AMB1051 

DO 1 I=1,KHM1 AMB1052 

MF=M AMB1053 

MHINV(I)=MINV0-DBLE(I-1)XDMINV AMB1054 

M=1.D0/MHINV(I) AMB1055 

DM=M-MF AMB1056 

M1=M-DM AMB1057 

M2=M-DM/2.D0 AMB1058 

M3=M AMB1059 

CALL MFUNC(M1,F1, ETALF1 , TETA1 ) AMB1060 

CALL MFUNC(M2,F2, ETALF2, TETA2) AMB1061 

CALL MFUNCCM3, F3,ETALF3,TETA3) AMB1062 

SUM=SUM+DMX(Fl+4.D0XF2+F3)/6 .DO AMB1063 

ETALF=ETALF3 AMB1064 

TETA=TETA3 AMB1065 

PSI=TETA+DASIN(1 . DO/M) AMB1066 

NORM=((3.DO-G)/4.DO)X(MXX2-1 . DO )XX0 . 75D0/ AMB1067 

1 (DSIN(PSI)X(l .D0+G1XMXX2)X*G14) AMB1068 

HM=SUMXNORM AMB1069 

HMV ( I ) =HM AMB1070 

GOREM=l . D0+G1XMXX2 AMB1071 

G0R=MX*2-1 . DO AMB1072 

DELT0B=0 . 5D0XDSQRT (GOR)X( 1 . DO/(MEXITXETALF) AMB1073 

1 +DSIN(TETA)/M)/DSIN( PSD+G15XHM/2 . DO AMB1074 

EPSIOB=DELTOB/DSQRT(GOR)-DSIN(TETA)/(MXDSIN(PSI)) AMB1075- 

KAPA0B=1 . DO AMB1076 

I F( DABS( PAI2-TETA ) .GT. 1 .D~6) AMB1077 

1KAPA0B=DTAN(TETA)XEPSI0B AMB1078 

LAMDOB=EPSIOB-DELTOBXGOREM/(GORXDSQRT(GOR)) AMB1079 

PRINT 11,I,M,HM,TETAXDEG,PSIXDEG AMB1080 
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11 F0RMAT(/1X,' I,M,HM,TETA,PSI=',I5,5D12.4) 

PRINT 12, DELTOB, EPSIOBXDEG, KAPAOB*DEG, LAMDOB*DEG 

12 FORMAT ( IX, * DELTOB, EPSI OB, KAPAO B, LAMD0B= ' , 5X, 5D12 . A) 

I CONTINUE 
MHINV(KHM)=0 . 

HMV(KHM) =1 . DO 
RETURN 

END 

SUBROUTINE MFUNCCM, F,ETALF,TETA) 

C SUBROUTINE NUMBER 21 

IMPLICIT REAL*8(A-H,0-Z,$) 

REAL*8 NU, NUFUNC, M, MEXIT, MD, MDD 

COMMON /GAMA/G, G1 ,G2,G3,GA,G5,G6,G7,G8,G9,G10,G11,G12,G13,G1A,G15 
1 G16,G17,G18,G19,G20 

COMMON /PAR/CO, ENO, EMI, D, SIGMA, TLIM,DRO, EL 0,Q0, TO, FACT, ALOGF, 

1 DPSIO,DTMAX,DETAO,ETALIM,XSI,XSF 

COMMON /GEOM/APF, PAI,PAI2,W,SH,CW,BETA,SBETA,CBETA,PSI1,SPSI1, 

1 CPSIl,PSIF,SPSIF,CPSIF,TPSIF,AK,SK,CK,AO,RF,XF,YF,ZF 

2 PHISOF,PHIF,SPHIF,CPHIF, DYMIN, RMIN,XS,DIST,XO,YO,ZO, 

3 DYO , DEG, PSIN, ST1 , CT1 , OMEGX, OMEGY, OMEGZ, XSVC 21 ) 

C 

QF(MDD)=1 . DO/DSQRTCMDDXX2-1 .DO) 

NUFUNC(MD)=-G5*DATAN(G5*QF(MD) )+DATAN(QF(MD) ) 

C 

MEXIT=EM1 

NU=NUFUNC(M) 

TETA= NUFUNC (MEXIT)+PA I 2-NU 
GOREM=1.DO+G1*M**2 
G0R=M**2-1 .DO 

F=(M**2)*(G0REM**G13)*DSIN(TETA)/G0R**1 .25D0 
GOREMl=l . DO+G1XMEXITXX2 
G0R1=MEXITXX2-1 .DO 

ETALF=( (GOREM/GOREM1 )x*GlA)x( (GORl/GOR)xxo . 25D0) 

RETURN 

END 

SUBROUTINE HINTERCM, H) 

C SUBROUTINE NUMBER 22 

IMPLICIT REALX8(A-H,0-Z,$) 

REAL*8 MINV,M, MEXIT, MHINV 

COMMON /GAMA/G, Gl, G2, G3, GA, G5, G6, G7, G8, G9, G1 0, Gil, G12, G13,G1 A, G15 
1 G16,G17,G18,G19,G20 

COMMON /PAR/CO, ENO, EMI, D, SIGMA, TLIM,DRO, EL 0,Q0, TO, FACT, ALOGF, 

1 DPSI 0 , DTMAX, DETAO, ETAL IM, XSI , XSF 

COMMON /G RP/ DM INV, MHINV (101), HMV( 101) 

COMMON /IGRP/KHM 
C COMPUTE H(M) BY INTERPOLATION 
MEXIT=EM1 
MINV=1 .DO/M 

I=KHM-IDINT(MINV/DMINV-1 .D-9)-l 
I F ( I .GE.l.AND. I .LT.KHM) GO TO 1 
PRINT 11,I,KHM,M,MEXIT 

II FORMAT (/IX, , I,KHM,M,MEXIT=*,2I5,2D1A.6/) 

CALL SOF( *2201 * ) 

1 CONTINUE 

F1=(MINV-MHINV(I+1))/DMINV 
F2=l . D0-F1 

IF(F1 . LT.-l . D-9 ) CALL SOFC2210') 

IF(F2.LT.-l.D-9) CALL S0FC2211*) 

H=F1XHMV(I)+F2XHMV(I+1) 

RETURN 

END 

SUBROUTINE MATCH( T, PSI 0, MAB , TETAAB) 

C SUBROUTINE NUMBER 23 

IMPLICIT REAL*8(A-H,0-Z,$) 

REAL *8 M, MOB, MEXIT, MAB, LAMDOB,KAPAOB 

COMMON /GAMA/G, Gl, G2, G3, GA, G5, G6 , G7 , G8, G9 , G1 0 , G1 1 , G12, G13 , G1 A, G15 
1 G16,G17,G18,G19,G20 

COMMON /PAR/CO, ENO, EMI, D, SIGMA, TLIM,DRO, EL 0,Q0, TO, FACT, ALOGF, 

1 DPSI 0, DTMAX, DETAO, ETAL IM, XSI, XSF 

COMMON /NPAR/NPHI,IPAR,NP,NR,NX,NXS,NS,NSPEC,NS1,NS2,NTAU0,NETA0, 
1 NAMB , NCASE, I CASE, I FAN 

COMMON /GEOM/APF, PAI , PAI2, H, SW, CH, BETA, SBETA, CBETA, PSI1 , SPSI 1 , 



AMB1081 
AMB1082 
AMB1083 
AMB108A 
AMB1085 
AMB1086 
AMB1087 
AMB1088 
AMB1089 
AMB1090 
AMB1091 
AMB1092 
AMB1093 
AMB109A 
AMB1095 
AMB1096 
AMB1097 
AMB1098 
AMB1099 
AMB1100 
AMB1101 
AMB1102 
AMB1103 
AMB110A 
AMB1105 
AMB1106 
AMB1107 
AMB1108 
AMB1109 
AMB1110 
AMB1111 
AMB1112 
AMB1113 
AMB111A 
AMB1115 
AMB1116 
AMB1117 
AMB1118 
AMB1119 
AMB1120 
AMB1121 
AMB1122 
AMB1123 
AMB112A 
AMB1125 
AMB1126 
AMB1127 
AMB 1128 
AMB1129 
AMB1130 
AMB 1131 
AMB 11 32 
AMB1133 
AMB113A 
AMB 1135 
AMB1136 
AMB 11 37 
AMB1138 
AMB1139 
AMB11A0 
AMB11A1 
AMB11A2 
AMB11A3 
AMB11AA 
AMB11A5 
AMB11A6 
AMB11A7 
AMB11A8 
AMB11A9 
AMB1150 
AMB1151 
AMB1152 
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1 CPSI1,PSIF,SPSIF,CPSIF,TPSIF,AK,SK,CK,A0,RF,XF,YF,ZF, 

2 PHISOF,PHIF,SPHIF,CPHIF,DYMIN,RMIN,XS,DIST,XO,YO,ZO, 

3 DYO,DEG,PSIN,ST1,CT1,OMEGX,OMEGY,OMEGZ,XSVC21) 

COMMON /POINT/XP,YP,XCOR,YCOR 

COMMON /GRP/DMINV ,MHINVC101), HMVC 1 01 ) 

COMMON /IGRP/KHM 
MEXIT=EM1 

GO TO C 1 01 , 102) , IFAN 

101 CONTINUE 

C FAN APPROXIMATED AS PLANAR 

MAB=DSQRTC1 . D0+G4/DTANC CPSI0-PSIF)/G5)**2) 

TETAAB=PSI0-DASINC1 .DO/MAB) 

GO TO 100 

102 CONTINUE 

C COMPUTE MAB FROM THE INVERSE PROBLEM SOLUTION 
COTAV=1.DO/DTAN(PSIO) 

EVY=YP*DLOG(YP/YCOR)/(YP-YCOR)-l . DO 

PSIN=PSI0 

DO 1 ITER=1,10 

PSI =PSIN 

M=DSQRT(1 .D0+G4/DTANCCPSI-PSIF)/G5)**2) 

M=DMAX1CM,MEXIT) 

CALL HINTERCM, HM) 

CALL MFUNCCM,F,ETALF,TETA) 

GOREM=l .D0+G1XMXX2 
G0R=M**2-1 . DO 

DELT0B=0.5D0XDSQRTCG0R)*C1. DO/ CMEXITxETAL F) 

1 +DSINCTETA)/M)/DSINCPSI)+G15*HM/2.D0 

EPSIOB=DELTOB/DSQRTCGOR)-DSINCTETA)/CM*DSINCPSI)) 

LAMDOB = EPSIOB-DELTOBXGOREM/C GORXDSQRTCGOR) ) 

COTN = COTAV+LAMDOB*EVY/DSINC PSI )X*2 
PSIN=PAI2-DATAN(C0TN) 

DPSI=PSIN-PSI 

IFCDABSCDPSD.LT. l.D-6) GO TO 11 

I CONTINUE 

PRINT 12,I,ITER,PSI,PSIN,DPSI,M,XP,YP,T 
12 FORMAT C/1X, ' I , ITER, PSI , PSIN, DPSI , M, XP, YP, T= •// 

1 1X,2I4,7D11 .3/) 

CALL SOFC ' 2301 ' ) 

II CONTINUE 

C USING M0B=M AS COMPUTED FROM THE INVERSE PROBLEM, FIND MAB. 

M0B = M 

CALL MFUNCCM, F,ETALF,TETA) 

PSI=TETA+DASIN(1 .DO/M) 

CALL HINTERCM, HM) 

GOREM=l . D0+G1XMXX2 
G0R=M**2-1 .DO 

DELT0B=0. 5DO*DSQRTCGOR)*C 1.D0/CMEXITXETALF) 

1 +DSINCTETA)/M)/DSINC PSI )+G15*HM/2 .DO 

FOB=CG7*GOREM)**G2/M 
FAB=F0B*C YP/YCOR) XXDELTOB 
CALL AREAFC FAB, MAB) 

EPSIOB=DELTOB/DSQRTCGOR)-DSINCTETA)/CM*DSINCPSI)) 

KAPA0B=1 .DO 

IFCDABSCPAI2-TETA) .GT.l.D-8) 

1KAPA0B=DTANCTETA)*EPSI0B 

COSTAB=DCOSCTETA)*CYP/YCOR)**C-KAPAOB) 

TETAAB = DACOSC COSTAB ) 

100 CONTINUE 
RETURN 
END 

SUBROUTINE AREAFCF,M) 

C SUBROUTINE NUMBER 24 

IMPLICIT REAL*8CA-H,0-Z,$) 

REALX8 MEXIT, MIN, M, MHINV 

COMMON /GAMA/G, G1 , G2,G3,G4,G5,G6,G7,G8,G9,G10,G11,G12,G13,G14,G15 
1 G16,G17,G18,G19,G20 

COMMON /PAR/CO, ENO, EMI, D, SIGMA, TLIM,DR0, EL 0,Q0, TO, FACT, ALOGF, 

1 DPSIO , DTMAX, DETAO , ETALIM, XSI , XSF 

COMMON /GRP/DMINV, MHINV Cl 01 ),HMVC 101) 

COMMON /IGRP/KHM 

C COMPUTE MACH NUMBER M FROM AREA RATIO FUNCTION F 



AMB1153 
AMB1154 
AMB1155 
AMB1156 
AMB1157 
AMB1158 
AMB1159 
AMB1160 
AMB1161 
AMB1162 
AMB1163 
AMB1164 
AMB1165 
AMB1166 
AMB1167 
AMB1168 
AMB1169 
AMB1170 
AMB1171 
AMB1172 
AMB1173 
AMB1174 
AMB1175 
AMB1176 
AMB1177 
AMB1178 
AMB1179 
AMB1180 
AMB1181 
AMB1182 
AMB1183 
AMB1184 
AMB1185 
AMB1186 
AMB1187 
AMB1188 
AMB1189 
AMB1190 
AMB1191 
AMB1192 
AMB1193 
AMB1194 
AMB1195 
AMB1196 
AMB1197 
AMB1198 
AMB1199 
AMB1200 
AMB1201 
AMB1202 
AMB1203 
AMB1204 
AMB1205 
AMB1206 
AMB1207 
AMB1208 
AMB1209 
AMB1210 
AMB1211 
AMB1212 
AMB1213 
AMB1214 
AMB1215 
AMB1216 
AMB1217 
. AMB1218 
AMB1219 
AMB1220 
AMB1221 
AMB1222 
AMB1223 
AMB1224 
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C F= 


(C2/(G+1))x(1+(G-1)xmx*2))xx((G+1)/(2x(G-1)))/M 


AMB1225 


C INITIAL GUESS IS MIN 


AMB1226 




MEXIT=EM1 


AMB1227 




E1=(FXMEXIT)XX(1 .D0/G2)/G7 


AMB1228 




E2=(E1-1.D0)/G1 


AMB1229 




E3=DMAX1(E2,MEXITXX2) 


AMB1230 




MIN=DSQRT(E3) 


AMB1231 




EMN=MIN 


AMB1232 




DO 1 1=1,100 


AMB1233 




EMO=EMN 


AMB1234 




GOREM=l . D0+G1 XEM0XX2 


AMB1235 




G0R=EM0XX2-1 .DO 


AMB1236 




F0=(G7 XGOREM) XXG2/ EMO 


AMB1237 




DF=FO-F 


AMB1238 


C 


PRINT 123,1, EMO, EMN, FO, F, DF,GOR,GOREM 


AMB1239 


C123 


FORMAT ( IX, ' I , EMO, EMN, F0,F,DF,G0R,G0REM=',I5,7D12.4) 


AMB1240 




DFDM=F1XG0R/(EM0XG0REM) 


AMB1241 




DMN=DF/DFDM 


AMB1242 




EMN=EMO-DMN 


AMB1243 




EPSEM= DABS ( DMN/ EMN) 


AMB1244 




IF(EPSEM.LT.l.D-lO) GO TO 11 


AMB1245 


1 


CONTINUE 


AMB1246 




CALL SOF( *2401 ' ) 


AMB1247 


11 


CONTINUE 


AMB1248 




M=EMN 


AMB1249 




RETURN 


AMB1250 




END 


AMB1251 
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